MODULE INIT      
 
! Allocate numerous arrays and perform miscellaneous initializations

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE MESH_POINTERS
USE GLOBAL_CONSTANTS
USE TRAN
USE MEMORY_FUNCTIONS, ONLY : CHKMEMERR 
USE DEVICE_VARIABLES
IMPLICIT NONE

PRIVATE
INTEGER IZERO
CHARACTER(255), PARAMETER :: initid='$Id: init.f90 6197 2010-05-11 13:20:05Z mcgratta $'
CHARACTER(255), PARAMETER :: initrev='$Revision: 6197 $'
CHARACTER(255), PARAMETER :: initdate='$Date: 2010-05-11 14:20:05 +0100 (Tue, 11 May 2010) $'

PUBLIC INITIALIZE_MESH_VARIABLES,INITIALIZE_GLOBAL_VARIABLES, OPEN_AND_CLOSE, INITIAL_NOISE, GET_REV_init

TYPE (MESH_TYPE), POINTER :: M
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (MATERIAL_TYPE), POINTER :: ML
 
 
CONTAINS
 
 
SUBROUTINE INITIALIZE_MESH_VARIABLES(NM)

USE PHYSICAL_FUNCTIONS, ONLY: GET_VISCOSITY,GET_SPECIFIC_GAS_CONSTANT,GET_SPECIFIC_HEAT
USE GEOMETRY_FUNCTIONS, ONLY: ASSIGN_PRESSURE_ZONE
USE RADCONS, ONLY: UIIDIM
USE CONTROL_VARIABLES
INTEGER :: N,I,J,K,II,JJ,KK,NN,IPTS,JPTS,KPTS,N_EDGES_DIM,IW,IC,IBC,IOR,IOPZ,IERR,IB,JB,KB,IPC,IPZ,ITMP
INTEGER, INTENT(IN) :: NM
REAL(EB) :: MU_N,YY_GET(1:N_SPECIES),VC,RTRM,CP_MF
INTEGER, POINTER :: IBP1, JBP1, KBP1,IBAR, JBAR, KBAR, NDWC, N_EDGES, NWC
REAL(EB),POINTER :: XS,XF,YS,YF,ZS,ZF
TYPE (INITIALIZATION_TYPE), POINTER :: IN
TYPE (P_ZONE_TYPE), POINTER :: PZ
TYPE (DEVICE_TYPE), POINTER :: DV
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC
 
IERR = 0
M => MESHES(NM)
IBP1 =>M%IBP1
JBP1 =>M%JBP1
KBP1 =>M%KBP1
IBAR =>M%IBAR
JBAR =>M%JBAR
KBAR =>M%KBAR
NDWC =>M%NDWC
NWC =>M%NWC
N_EDGES=>M%N_EDGES
XS=>M%XS 
YS=>M%YS 
ZS=>M%ZS
XF=>M%XF 
YF=>M%YF 
ZF=>M%ZF
ALLOCATE(  M%RHO(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHO',IZERO)
ALLOCATE( M%RHOS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHOS',IZERO)
ALLOCATE(  M%TMP(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','TMP',IZERO)
ALLOCATE( M%FRHO(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FRHO',IZERO)

ALLOCATE( M%U(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','U',IZERO)
ALLOCATE( M%V(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','V',IZERO)
ALLOCATE( M%W(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','W',IZERO)
ALLOCATE(M%US(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','US',IZERO)
ALLOCATE(M%VS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','VS',IZERO)
ALLOCATE(M%WS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WS',IZERO)
ALLOCATE(M%FVX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVX',IZERO)
ALLOCATE(M%FVY(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVY',IZERO)
ALLOCATE(M%FVZ(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVZ',IZERO)
ALLOCATE(   M%H(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','H',IZERO)
ALLOCATE(   M%HS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','HS',IZERO)
ALLOCATE(   M%KRES(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','KRES',IZERO)
ALLOCATE(M%DDDT(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','DDDT',IZERO)
ALLOCATE(M%D(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','D',IZERO)
ALLOCATE(M%DS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','DS',IZERO)
ALLOCATE(M%MU(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','MU',IZERO)
IF (STORE_MU_DNS .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%MU_DNS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','MU_DNS',IZERO)
ENDIF
IF (.NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%Q(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','Q',IZERO)
ENDIF

IF (DYNSMAG .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%C_DYNSMAG(IBAR,JBAR,KBAR),STAT=IZERO)
   CALL ChkMemErr('INIT','C_DYNSMAG',IZERO)
ENDIF
IF (CHECK_KINETIC_ENERGY .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ! measure of turbulence resolution
   ALLOCATE(M%MTR(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','MTR',IZERO)
   M%MTR=0._EB
   ! measure of scalar resolution
   ALLOCATE(M%MSR(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','MSR',IZERO)
   M%MSR=0._EB
   ! wavelet error measure
   ALLOCATE(M%WEM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','WEM',IZERO)
   M%WEM=0._EB
ENDIF

! Background pressure, temperature, density as a function of height (Z coordinate)

ALLOCATE(  M%PBAR(0:KBP1,0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','PBAR',IZERO)
ALLOCATE(  M%PBAR_S(0:KBP1,0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','PBAR_S',IZERO)
ALLOCATE(  M%R_PBAR(0:KBP1,0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','R_PBAR',IZERO)
ALLOCATE(  M%D_PBAR_DT(N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','D_PBAR_DT',IZERO)
ALLOCATE(  M%D_PBAR_S_DT(N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','D_PBAR_S_DT',IZERO)
ALLOCATE(M%P_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','P_0',IZERO)
ALLOCATE(M%TMP_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','TMP_0',IZERO)
ALLOCATE(M%RHO_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHO_0',IZERO)

! Leaks

ALLOCATE(  M%U_LEAK(0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','U_LEAK',IZERO)
M%U_LEAK = 0._EB

! Allocate species arrays
 
IF (N_SPECIES>0 .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE( M%YY(0:IBP1,0:JBP1,0:KBP1,N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','YY',IZERO)
   M%YY = 0._EB
   ALLOCATE(M%YYS(0:IBP1,0:JBP1,0:KBP1,N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','YYS',IZERO)
   M%YYS = 0._EB   
   ALLOCATE(M%DEL_RHO_D_DEL_Y(0:IBP1,0:JBP1,0:KBP1,N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','DEL_RHO_D_DEL_Y',IZERO)
   M%DEL_RHO_D_DEL_Y = 0._EB   
   ALLOCATE(M%RSUM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','RSUM',IZERO)
   M%RSUM = RSUM0
ENDIF

! Allocate water mass arrays if sprinklers are present
 
IF (DROPLET_FILE) PARTICLE_TAG = NM

IF (N_EVAP_INDICES>0 .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%QR_W(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','QR_W',IZERO) 
   M%QR_W = 0._EB
   ALLOCATE(M%AVG_DROP_DEN(0:IBP1,0:JBP1,0:KBP1,N_EVAP_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_DEN',IZERO) 
   M%AVG_DROP_DEN=0._EB
   ALLOCATE(M%AVG_DROP_AREA(0:IBP1,0:JBP1,0:KBP1,N_EVAP_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_AREA',IZERO) 
   M%AVG_DROP_AREA=0._EB
   ALLOCATE(M%AVG_DROP_TMP(0:IBP1,0:JBP1,0:KBP1,N_EVAP_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_TMP',IZERO) 
   M%AVG_DROP_TMP=TMPM
   ALLOCATE(M%AVG_DROP_RAD(0:IBP1,0:JBP1,0:KBP1,N_EVAP_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_RAD',IZERO) 
   M%AVG_DROP_RAD=0._EB
   ALLOCATE(M%D_LAGRANGIAN(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','D_LAGRANGIAN',IZERO) 
   M%D_LAGRANGIAN = 0._EB
ENDIF
 
! If radiation absorption desired allocate arrays
 
IF (.NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%QR(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','QR',IZERO)
   M%QR = 0._EB
   ALLOCATE(M%KAPPA(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','KAPPA',IZERO) 
   M%KAPPA = KAPPA0
   ALLOCATE(M%UII(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','UII',IZERO)
   M%UII = 0._EB
ELSE
   ALLOCATE(M%QR(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','QR',IZERO)
   M%QR = 0._EB
ENDIF

! Work arrays

ALLOCATE(M%WORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK1',IZERO)
ALLOCATE(M%WORK2(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK2',IZERO)
ALLOCATE(M%WORK3(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK3',IZERO)
ALLOCATE(M%WORK4(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK4',IZERO)
ALLOCATE(M%WORK5(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK5',IZERO)
ALLOCATE(M%WORK6(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK6',IZERO)
ALLOCATE(M%WORK7(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WORK7',IZERO)

IF (FLUX_LIMITER>=0 .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%SCALAR_SAVE1(0:IBP1,0:JBP1,0:KBP1,0:N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','SCALAR_SAVE1',IZERO)
   ALLOCATE(M%SCALAR_SAVE2(0:IBP1,0:JBP1,0:KBP1,0:N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','SCALAR_SAVE2',IZERO)
   ALLOCATE(M%SCALAR_SAVE3(0:IBP1,0:JBP1,0:KBP1,0:N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','SCALAR_SAVE3',IZERO)
   IF (N_SPECIES>0) THEN
      ALLOCATE(M%SCALAR_SAVE4(0:IBP1,0:JBP1,0:KBP1,1:N_SPECIES),STAT=IZERO)
      CALL ChkMemErr('INIT','SCALAR_SAVE4',IZERO)
      ALLOCATE(M%SCALAR_SAVE5(0:IBP1,0:JBP1,0:KBP1,1:N_SPECIES),STAT=IZERO)
      CALL ChkMemErr('INIT','SCALAR_SAVE5',IZERO)
   ENDIF
ENDIF

IF (IMMERSED_BOUNDARY_METHOD==2) THEN
   ALLOCATE(M%IBM_SAVE1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','IBM_SAVE1',IZERO)
   ALLOCATE(M%IBM_SAVE2(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','IBM_SAVE2',IZERO)
   ALLOCATE(M%IBM_SAVE3(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','IBM_SAVE3',IZERO)
   ALLOCATE(M%IBM_SAVE4(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','IBM_SAVE4',IZERO)
   ALLOCATE(M%IBM_SAVE5(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','IBM_SAVE5',IZERO)
   ALLOCATE(M%IBM_SAVE6(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','IBM_SAVE6',IZERO)
ENDIF

! Boundary file patch counter
 
ALLOCATE(M%INC(-3:3,0:M%N_OBST),STAT=IZERO)
CALL ChkMemErr('INIT','INC',IZERO)

! Initialize background pressure, temperature and density

M%D_PBAR_DT   = 0._EB
M%D_PBAR_S_DT = 0._EB

IF (STRATIFICATION .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   DO K=0,M%KBP1
      M%TMP_0(K) = TMPA + LAPSE_RATE*(M%ZC(K)-GROUND_LEVEL)
      IF (LAPSE_RATE /= 0._EB) THEN
         M%P_0(K) = P_INF*(M%TMP_0(K)/M%TMP_0(0))**(GVEC(3)/RSUM0/LAPSE_RATE)
      ELSE
         M%P_0(K) = P_INF*EXP(GVEC(3)*(M%ZC(K)-GROUND_LEVEL)/(RSUM0*TMPA))
      ENDIF
   ENDDO
ELSE
   M%TMP_0(:) = TMPA
   M%P_0(:)   = P_INF
ENDIF
DO K=0,M%KBP1
   M%PBAR(K,:)   = M%P_0(K)
   M%PBAR_S(K,:) = M%P_0(K)
   M%RHO_0(K)    = M%P_0(K)/(M%TMP_0(K)*RSUM0)
ENDDO

! Initialize various time step variables
 
M%DT_PREV = M%DT
M%DT_NEXT = M%DT
M%DT_INIT = M%DT

! Initialize major arrays

DO K=0,M%KBP1
   M%RHO(:,:,K) = M%RHO_0(K)
   M%TMP(:,:,K) = M%TMP_0(K)
ENDDO
M%RHOS    = M%RHO
IF (.NOT.EVACUATION_ONLY(NM)) M%FRHO    = 0._EB
M%U       = U0
M%V       = V0
M%W       = W0
M%US      = U0
M%VS      = V0
M%WS      = W0
M%FVX   = 0._EB
M%FVY   = 0._EB
M%FVZ   = 0._EB
M%H     = H0
M%HS    = H0
M%KRES  = 0._EB

M%DDDT  = 0._EB
M%D     = 0._EB
M%DS    = 0._EB
IF (.NOT.EVACUATION_ONLY(NM)) THEN
   M%Q     = 0._EB
ENDIF
IF (EVACUATION_ONLY(NM)) THEN
   M%U       = 0._EB
   M%V       = 0._EB
   M%W       = 0._EB
   M%US      = 0._EB
   M%VS      = 0._EB
   M%WS      = 0._EB
   M%H       = 0._EB
ENDIF
IF (N_SPECIES > 0 .AND. .NOT.EVACUATION_ONLY(NM)) M%DEL_RHO_D_DEL_Y = 0._EB
! Viscosity
IF (N_SPECIES > 0) THEN
    YY_GET(:) = SPECIES(1:N_SPECIES)%YY0
   CALL GET_VISCOSITY(YY_GET,MU_N,NINT(TMPA))
ELSE
   MU_N = Y2MU_C(NINT(TMPA))*SPECIES(0)%MW
ENDIF

M%MU = MU_N
IF (STORE_MU_DNS) M%MU_DNS = MU_N

IF (DYNSMAG .AND. .NOT.EVACUATION_ONLY(NM)) M%C_DYNSMAG = 0._EB
 
! Initialize mass fraction arrays

IF (N_SPECIES > 0 .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   DO N=1,N_SPECIES
      M%YY(:,:,:,N)  = SPECIES(N)%YY0
      M%YYS(:,:,:,N) = SPECIES(N)%YY0
   ENDDO
ENDIF

! Initialize pressure ZONEs
 
IF (.NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%PRESSURE_ZONE(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','PRESSURE_ZONE',IZERO)
   M%PRESSURE_ZONE = 0
ENDIF
ZONE_LOOP: DO N=1,N_ZONE
   IF (EVACUATION_ONLY(NM)) CYCLE ZONE_LOOP
   PZ => P_ZONE(N)
   DO K=0,KBP1
      DO J=0,JBP1
         DO I=0,IBP1
            IF (M%PRESSURE_ZONE(I,J,K)==N) CYCLE
            IF (M%XC(I) > PZ%X1 .AND. M%XC(I) < PZ%X2 .AND. &
                M%YC(J) > PZ%Y1 .AND. M%YC(J) < PZ%Y2 .AND. &
                M%ZC(K) > PZ%Z1 .AND. M%ZC(K) < PZ%Z2) THEN 
                M%PRESSURE_ZONE(I,J,K) = N
                DO IOPZ=0,N_ZONE
                   IF (PZ%LEAK_AREA(IOPZ) > 0._EB) ACTUAL_LEAK_AREA(N,IOPZ) = PZ%LEAK_AREA(IOPZ)
                   IF (PZ%LEAK_AREA(IOPZ) > 0._EB) ACTUAL_LEAK_AREA(IOPZ,N) = PZ%LEAK_AREA(IOPZ)
                ENDDO
                IF (.NOT.M%SOLID(M%CELL_INDEX(I,J,K))) CALL ASSIGN_PRESSURE_ZONE(NM,M%XC(I),M%YC(J),M%ZC(K),N)
            ENDIF
         ENDDO
      ENDDO
   ENDDO
ENDDO ZONE_LOOP


! Over-ride default ambient conditions with user-prescribed INITializations
 
DO N=1,N_INIT
   IF (EVACUATION_ONLY(NM)) CYCLE
   IN => INITIALIZATION(N)
   DO K=0,KBP1
      DO J=0,JBP1
         DO I=0,IBP1
            IF (M%XC(I) > IN%X1 .AND. M%XC(I) < IN%X2 .AND. &
                M%YC(J) > IN%Y1 .AND. M%YC(J) < IN%Y2 .AND. &
                M%ZC(K) > IN%Z1 .AND. M%ZC(K) < IN%Z2) THEN
               M%TMP(I,J,K)            = IN%TEMPERATURE
               M%RHO(I,J,K)            = IN%DENSITY
               IF (N_SPECIES>0) M%YY(I,J,K,1:N_SPECIES) = IN%MASS_FRACTION(1:N_SPECIES)
               IF (IN%ADJUST_DENSITY)     M%RHO(I,J,K) = M%RHO(I,J,K)*M%P_0(K)/P_INF
               IF (IN%ADJUST_TEMPERATURE) M%TMP(I,J,K) = M%TMP(I,J,K)*M%P_0(K)/P_INF
            ENDIF
         ENDDO
      ENDDO
   ENDDO
ENDDO

M%RHOS    = M%RHO

! Compute molecular weight term RSUM=R0*SUM(Y_i/M_i)

IF (N_SPECIES>0 .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            YY_GET(:) = M%YY(I,J,K,:)
            CALL GET_SPECIFIC_GAS_CONSTANT(YY_GET,M%RSUM(I,J,K))
         ENDDO
      ENDDO
   ENDDO
ENDIF

! Allocate and Initialize Mesh-Dependent Radiation Arrays

M%QR    = 0._EB
IF (.NOT.EVACUATION_ONLY(NM)) THEN
   M%KAPPA = KAPPA0
   M%UII   = 4._EB*SIGMA*TMPA4
ENDIF
M%ANGLE_INC_COUNTER = 0
M%RAD_CALL_COUNTER  = 0
IF (RADIATION .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%UIID(0:M%IBP1,0:M%JBP1,0:M%KBP1,1:UIIDIM),STAT=IZERO)
   CALL ChkMemErr('INIT','UIID',IZERO)
   M%UIID = 0.
ENDIF

! General work arrays

M%WORK1 = 0._EB
M%WORK2 = 0._EB
M%WORK3 = 0._EB
M%WORK4 = 0._EB
M%WORK5 = 0._EB
M%WORK6 = 0._EB
IF (.NOT.EVACUATION_ONLY(NM)) M%WORK7 = 0._EB

! Scalar

IF (FLUX_LIMITER>=0) THEN
   M%SCALAR_SAVE1 = 0._EB
   M%SCALAR_SAVE2 = 0._EB
   M%SCALAR_SAVE3 = 0._EB
   IF (N_SPECIES>0) THEN
      M%SCALAR_SAVE4 = 0._EB
      M%SCALAR_SAVE5 = 0._EB
   ENDIF
ENDIF

! Immersed Boundary Method

IF (IMMERSED_BOUNDARY_METHOD==2) THEN
   M%IBM_SAVE1 = 0._EB
   M%IBM_SAVE2 = 0._EB
   M%IBM_SAVE3 = 0._EB
   M%IBM_SAVE4 = 0._EB
   M%IBM_SAVE5 = 0._EB
   M%IBM_SAVE6 = 0._EB
ENDIF

! Designate each boundary cell with a reference number for wall BC's
 
NWC  = 0
NDWC = 0
 
! Determine the number of wall cells to allocate
 
OBST_LOOP_1: DO N=1,M%N_OBST
   OB=>M%OBSTRUCTION(N)
   IF (OB%CONSUMABLE .AND. .NOT.EVACUATION_ONLY(NM)) THEN
      IB   = OB%I2-OB%I1
      JB   = OB%J2-OB%J1
      KB   = OB%K2-OB%K1
      NDWC = NDWC + 2*(MAX(1,IB)*JB*KB + MAX(1,JB)*IB*KB + MAX(1,KB)*IB*JB)
   ELSE
      DO K=OB%K1+1,OB%K2
         DO J=OB%J1+1,OB%J2
            IC = M%CELL_INDEX(OB%I1  ,J,K)
            IF (.NOT.M%SOLID(IC) .OR. M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) NDWC = NDWC + 1
            IC = M%CELL_INDEX(OB%I2+1,J,K)
            IF (.NOT.M%SOLID(IC) .OR. M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) NDWC = NDWC + 1
         ENDDO 
      ENDDO
      DO K=OB%K1+1,OB%K2
         DO I=OB%I1+1,OB%I2
            IC = M%CELL_INDEX(I,OB%J1  ,K)
            IF (.NOT.M%SOLID(IC) .OR. M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) NDWC = NDWC + 1
            IC = M%CELL_INDEX(I,OB%J2+1,K)
            IF (.NOT.M%SOLID(IC) .OR. M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) NDWC = NDWC + 1
         ENDDO 
      ENDDO
      DO J=OB%J1+1,OB%J2
         DO I=OB%I1+1,OB%I2
            IC = M%CELL_INDEX(I,J,OB%K1  )
            IF (.NOT.M%SOLID(IC) .OR. M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) NDWC = NDWC + 1
            IC = M%CELL_INDEX(I,J,OB%K2+1)
            IF (.NOT.M%SOLID(IC) .OR. M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) NDWC = NDWC + 1
         ENDDO 
      ENDDO
   ENDIF
ENDDO OBST_LOOP_1

! Add wall cells for VIRTUAL devices

M%NVWC = 0
DO N=1,N_DEVC
   DV => DEVICE(N)
   IF (DV%MESH/=NM) CYCLE
   IF (EVACUATION_ONLY(NM)) CYCLE
   IF (DV%QUANTITY/='CABLE TEMPERATURE') CYCLE
   M%NVWC = M%NVWC + 1
ENDDO

DO N=1,N_PART
   IF (EVACUATION_ONLY(NM)) CYCLE
   PC => PARTICLE_CLASS(N)
   IF (PC%SURF_INDEX>0) THEN
      DO IB=1,N_INIT
         IN => INITIALIZATION(IB)
         IF (IN%PART_INDEX==N) M%NVWC = M%NVWC + IN%NUMBER_INITIAL_DROPLETS
      ENDDO
   ENDIF
ENDDO

! Allocate arrays indexed by wall cells (IW)

NDWC = NDWC + M%NEWC + M%NVWC

ALLOCATE(M%WALL(1:NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','WALL',IZERO)

ALLOCATE(M%RHO_F(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','RHO_F',IZERO) 
M%RHO_F = RHOA
ALLOCATE(M%E_WALL(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','E_WALL',IZERO)
M%E_WALL = 1._EB
ALLOCATE(M%U_TAU(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','U_TAU',IZERO)
M%U_TAU = 0._EB

IF (.NOT.EVACUATION_ONLY(NM)) THEN 

   ALLOCATE(M%TMP_F(NDWC),STAT=IZERO)
   CALL ChkMemErr('INIT','TMP_F',IZERO) 
   M%TMP_F = TMPA
   ALLOCATE(M%TMP_B(NDWC),STAT=IZERO)
   CALL ChkMemErr('INIT','TMP_B',IZERO) 
   M%TMP_B = TMPA
   ALLOCATE(M%YY_F(NDWC,N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','YY_F',IZERO)  
   DO IW=1,NDWC
      M%YY_F(IW,1:N_SPECIES) = SPECIES(1:N_SPECIES)%YY0
   ENDDO
   ALLOCATE(M%QRADIN(NDWC),STAT=IZERO)
   CALL ChkMemErr('INIT','QRADIN',IZERO) 
   M%QRADIN = SIGMA*TMPA4
   ALLOCATE(M%QRADOUT(NDWC),STAT=IZERO)
   CALL ChkMemErr('INIT','QRADOUT',IZERO) 
   M%QRADOUT = SIGMA*TMPA4
   ALLOCATE(M%QCONF(NDWC),STAT=IZERO)
   CALL ChkMemErr('INIT','QCONF',IZERO) 
   M%QCONF = 0._EB
   ALLOCATE(M%HEAT_TRANS_COEF(NDWC),STAT=IZERO)
   CALL ChkMemErr('INIT','HEAT_TRANS_COEF',IZERO) 
   M%HEAT_TRANS_COEF = 0._EB
   
   ALLOCATE(M%D_CORR(M%NEWC),STAT=IZERO)
   CALL ChkMemErr('INIT','D_CORR',IZERO) 
   ALLOCATE(M%DS_CORR(M%NEWC),STAT=IZERO)
   CALL ChkMemErr('INIT','DS_CORR',IZERO) 
   M%D_CORR = 0._EB
   M%DS_CORR = 0._EB
   ALLOCATE(M%UVW_SAVE(M%NEWC),STAT=IZERO)
   CALL ChkMemErr('INIT','UVW_SAVE',IZERO) 
   M%UVW_SAVE = 0._EB
   
   IF (OLD_PRESSURE_SCHEME .AND. ITERATE_PRESSURE) THEN
      ALLOCATE(M%DH(M%NEWC),STAT=IZERO)
      CALL ChkMemErr('INIT','DH',IZERO)
      M%DH=0._EB
      ALLOCATE(M%DHS(M%NEWC),STAT=IZERO)
      CALL ChkMemErr('INIT','DHS',IZERO)
      M%DHS=0._EB
   ENDIF

ENDIF

IF (FLUX_LIMITER>=0) THEN
   ! diffusive flux at the wall
   ALLOCATE(M%FW(NDWC,N_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','FW',IZERO)  
   DO IW=1,NDWC
      M%FW(IW,1:N_SPECIES) = 0._EB
   ENDDO
ENDIF

ALLOCATE(M%XW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','XW',IZERO)
ALLOCATE(M%YW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','YW',IZERO)
ALLOCATE(M%ZW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','ZW',IZERO)
ALLOCATE(M%TW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','TW',IZERO)  
M%TW = T_BEGIN

IF (.NOT.EVACUATION_ONLY(NM)) THEN 

ALLOCATE(M%EW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','EW',IZERO)  
M%EW = 0._EB
ALLOCATE(M%KW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','KW',IZERO)  
M%KW = 1._EB
ALLOCATE(M%RHODW(NDWC,N_SPECIES),STAT=IZERO)
CALL ChkMemErr('INIT','RHODW',IZERO) 
M%RHODW = 1._EB
ALLOCATE(M%AREA_ADJUST(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','AREA_ADJUST',IZERO)  
M%AREA_ADJUST = 1._EB
ALLOCATE(M%MASSFLUX(NDWC,0:N_SPECIES),STAT=IZERO)
CALL ChkMemErr('INIT','MASSFLUX',IZERO)  
M%MASSFLUX = 0._EB
ALLOCATE(M%MASSFLUX_ACTUAL(NDWC,0:N_SPECIES),STAT=IZERO)
CALL ChkMemErr('INIT','MASSFLUX_ACTUAL',IZERO)  
M%MASSFLUX_ACTUAL = 0._EB
ALLOCATE(M%NPPCW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','NPPCW',IZERO) 
M%NPPCW = 1
ALLOCATE(M%PRESSURE_ZONE_WALL(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','PRESSURE_ZONE_WALL',IZERO) 
M%PRESSURE_ZONE_WALL = 0
ALLOCATE(M%WALL_INDEX_BACK(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_INDEX_BACK',IZERO) 
M%WALL_INDEX_BACK = 0
ALLOCATE(M%AW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','AW',IZERO)  
M%AW = 0._EB
ALLOCATE(M%RAW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','RAW',IZERO) 
M%RAW = 0._EB

ENDIF

ALLOCATE(M%RDN(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','RDN',IZERO)
M%RDN = 1._EB
ALLOCATE(M%UW0(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','UW0',IZERO)    
M%UW0 = 0._EB
ALLOCATE(M%UW(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','UW',IZERO)    
M%UW = 0._EB
ALLOCATE(M%UWS(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','UWS',IZERO)   
M%UWS = 0._EB
ALLOCATE(M%OBST_INDEX_W(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','OBST_INDEX_W',IZERO) 
M%OBST_INDEX_W = 0
ALLOCATE(M%VENT_INDEX(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','VENT_INDEX',IZERO) 
M%VENT_INDEX = 0
ALLOCATE(M%DUWDT(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','DUWDT',IZERO) 
M%DUWDT = 0._EB
ALLOCATE(M%PRESSURE_BC_INDEX(M%NEWC),STAT=IZERO)
CALL ChkMemErr('INIT','PRESSURE_BC_INDEX',IZERO) 
M%PRESSURE_BC_INDEX = NEUMANN
ALLOCATE(M%IJKW(15,NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','IJKW',IZERO)  
M%IJKW = 0
ALLOCATE(M%IBC_ORIG(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','IBC_ORIG',IZERO)  
M%IBC_ORIG = 0
ALLOCATE(M%BOUNDARY_TYPE(0:NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','BOUNDARY_TYPE',IZERO) 
M%BOUNDARY_TYPE = NULL_BOUNDARY
ALLOCATE(M%WALL_INDEX(0:M%CELL_COUNT,-3:3),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_INDEX',IZERO) 
M%WALL_INDEX = 0
ALLOCATE(M%EDGE_INDEX(0:M%CELL_COUNT,1:12),STAT=IZERO)
CALL ChkMemErr('INIT','EDGE_INDEX',IZERO) 
M%EDGE_INDEX = 0
ALLOCATE(M%UVW_GHOST(0:M%CELL_COUNT,3),STAT=IZERO)
CALL ChkMemErr('INIT','UVW_GHOST',IZERO) 
M%UVW_GHOST = 0

! Surface soot array

IF (SOOT_DEPOSITION .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%AWMSOOT(NDWC),STAT=IZERO)
   CALL ChkMemErr('INIT','AWMSOOT',IZERO) 
   M%AWMSOOT = 0._EB
ENDIF

! Surface water arrays
 
IF (ACCUMULATE_WATER .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%AWMPUA(NDWC,N_EVAP_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AWMPUA',IZERO) 
   M%AWMPUA = 0._EB
ENDIF

IF (.NOT.EVACUATION_ONLY(NM)) THEN

ALLOCATE(M%WMPUA(NDWC,MAX(1,N_EVAP_INDICES)),STAT=IZERO)
CALL ChkMemErr('INIT','WMPUA',IZERO) 
M%WMPUA = 0._EB
ALLOCATE(M%WCPUA(NDWC,MAX(1,N_EVAP_INDICES)),STAT=IZERO)
CALL ChkMemErr('INIT','WCPUA',IZERO) 
M%WCPUA = 0._EB

ENDIF
 
! Surface work arrays
 
ALLOCATE(M%WALL_WORK1(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_WORK1',IZERO) 
ALLOCATE(M%WALL_WORK2(NDWC),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_WORK2',IZERO) 

! Vegetation surface drag
ALLOCATE(M%VEG_DRAG(0:IBP1,0:JBP1),STAT=IZERO)
CALL ChkMemErr('INIT','VEG_DRAG',IZERO) 
M%VEG_DRAG = 0._EB
 
! Set up boundary arrays for external boundaries of the current mesh

DO K=1,KBAR
   DO J=1,JBAR
      I   = 0
      IBC = DEFAULT_SURF_INDEX
      IOR = 1
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC,IERR)
      IF (IERR>0) RETURN
   ENDDO
ENDDO
DO K=1,KBAR
   DO J=1,JBAR
      I   = IBP1
      IBC = DEFAULT_SURF_INDEX
      IOR = -1
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC,IERR)
      IF (IERR>0) RETURN
   ENDDO
ENDDO
 
DO K=1,KBAR
   DO I=1,IBAR
      J   = 0
      IBC = DEFAULT_SURF_INDEX
      IOR = 2
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC,IERR)
      IF (IERR>0) RETURN
   ENDDO
ENDDO
DO K=1,KBAR
   DO I=1,IBAR
      J   = JBP1
      IBC = DEFAULT_SURF_INDEX
      IOR = -2
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC,IERR)
      IF (IERR>0) RETURN
   ENDDO
ENDDO
 
IF (.NOT.EVACUATION_ONLY(NM)) THEN
DO J=1,JBAR
   DO I=1,IBAR
      K   = 0
      IBC = DEFAULT_SURF_INDEX
      IOR = 3
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC,IERR)
      IF (IERR>0) RETURN
   ENDDO
ENDDO
DO J=1,JBAR
   DO I=1,IBAR
      K   = KBP1
      IBC = DEFAULT_SURF_INDEX
      IOR = -3
      NWC = NWC + 1
      IW  = NWC
      CALL INIT_WALL_CELL(NM,I,J,K,0,IW,IOR,IBC,IERR)
      IF (IERR>0) RETURN
   ENDDO
ENDDO
ENDIF

! Go through all obstructions and decide which cell faces ought to be given a wall cell index and initialized
 
OBST_LOOP_2: DO N=1,M%N_OBST
   OB=>M%OBSTRUCTION(N)
 
   DO K=OB%K1+1,OB%K2
      DO J=OB%J1+1,OB%J2
         I = OB%I1+1 
         IF (I==1) CYCLE   ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I-1,J,K)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = -1 
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC,IERR)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO
 
   DO K=OB%K1+1,OB%K2
      DO J=OB%J1+1,OB%J2
         I = OB%I2
         IF (I==M%IBAR) CYCLE  ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I+1,J,K)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = 1
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC,IERR)
         IF (IERR>0) RETURN
      ENDDO 
   ENDDO
 
   DO K=OB%K1+1,OB%K2
      DO I=OB%I1+1,OB%I2
         J = OB%J1+1
         IF (J==1) CYCLE   ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I,J-1,K)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = -2
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC,IERR)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO   
 
   DO K=OB%K1+1,OB%K2
      DO I=OB%I1+1,OB%I2
         J = OB%J2
         IF (J==M%JBAR) CYCLE  ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I,J+1,K)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = 2
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC,IERR)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO   
 
   DO J=OB%J1+1,OB%J2
      DO I=OB%I1+1,OB%I2
         K = OB%K1+1
         IF (K==1) CYCLE   ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I,J,K-1)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = -3
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC,IERR)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO   
 
   DO J=OB%J1+1,OB%J2
      DO I=OB%I1+1,OB%I2
         K = OB%K2
         IF (K==M%KBAR) CYCLE  ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IC = M%CELL_INDEX(I,J,K+1)
         IF (M%SOLID(IC) .AND. .NOT.M%OBSTRUCTION(M%OBST_INDEX_C(IC))%REMOVABLE) CYCLE   ! Permanently covered face
         IOR = 3 
         IBC = OB%IBC(IOR)
         IW  = M%WALL_INDEX(IC,-IOR)
         IF (IW==0) THEN
            NWC = NWC + 1
            IW  = NWC
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,N,IW,IOR,IBC,IERR)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO   
 
ENDDO OBST_LOOP_2

!Initialize PSUM for zone cases
IF (N_ZONE > 0) THEN
   N_ZONE_LOOP: DO IPZ = 1,N_ZONE
      IF (EVACUATION_ONLY(NM)) EXIT N_ZONE_LOOP
      PSUM(IPZ,NM) = 0._EB
      DO K=1,M%KBAR
         DO J=1,M%JBAR
            DO I=1,M%IBAR
               IF (M%PRESSURE_ZONE(I,J,K) /= IPZ) CYCLE
               IF (M%SOLID(M%CELL_INDEX(I,J,K)))    CYCLE
               VC   = M%DX(I)*M%RC(I)*M%DY(J)*M%DZ(K)
               ITMP = MIN(5000,NINT(M%TMP(I,J,K)))
               IF (N_SPECIES==0) THEN
                  RTRM = RSUM0/(Y2CP_C(ITMP)*M%PBAR(K,IPZ))
               ELSE
                  YY_GET(:) = M%YY(I,J,K,:)
                  CALL GET_SPECIFIC_HEAT(YY_GET,CP_MF,ITMP)
                  RTRM = M%RSUM(I,J,K)/(CP_MF*M%PBAR(K,IPZ))
               ENDIF
               PSUM(IPZ,NM) = PSUM(IPZ,NM) + VC*(1._EB/M%PBAR(K,IPZ)-RTRM)
            ENDDO
         ENDDO
      ENDDO
   ENDDO N_ZONE_LOOP
ENDIF

! Set up wall cell arrays for VIRTUAL boundaries

IW = M%NWC

DEVICE_LOOP: DO N=1,N_DEVC
   IF (EVACUATION_ONLY(NM)) CYCLE DEVICE_LOOP
   DV => DEVICE(N)
   IF (DV%MESH/=NM) CYCLE DEVICE_LOOP
   IF (DV%QUANTITY/='CABLE TEMPERATURE') CYCLE DEVICE_LOOP
   IW = IW + 1
   DV%VIRTUAL_WALL_INDEX = IW 
   I = DV%I
   J = DV%J
   K = DV%K
   IBC = DV%SURF_INDEX
   CALL INIT_WALL_CELL(NM,I,J,K,0,IW,0,IBC,IERR)
ENDDO DEVICE_LOOP

INITIALIZATION_LOOP: DO N=1,N_INIT
   IF (EVACUATION_ONLY(NM)) CYCLE INITIALIZATION_LOOP
   IN => INITIALIZATION(N)
   IPC = IN%PART_INDEX
   IF (IPC<1) CYCLE INITIALIZATION_LOOP
   PC =>PARTICLE_CLASS(IPC)
   IF (PC%SURF_INDEX>0) THEN
      IN%WALL_INDEX_START = IW + 1
      DO NN=1,IN%NUMBER_INITIAL_DROPLETS
         IW = IW + 1
         I = 1
         J = 1
         K = 1
         IBC = PC%SURF_INDEX
         CALL INIT_WALL_CELL(NM,I,J,K,0,IW,0,IBC,IERR)
      ENDDO
   ENDIF
ENDDO INITIALIZATION_LOOP

! Determine back wall index for exposed surfaces

DO IW=M%NEWC+1,M%NWC  
   IF (EVACUATION_ONLY(NM)) CYCLE
   ! Only assign WALL_INDEX_BACK to wall cells that are not attached to the exterior boundary of the computational domain
   SF=>SURFACE(M%IJKW(5,IW))
   IF (SF%BACKING==EXPOSED) THEN
      II = M%IJKW(1,IW)
      JJ = M%IJKW(2,IW)
      KK = M%IJKW(3,IW)
      IC = M%CELL_INDEX(II,JJ,KK)
      IOR = M%IJKW(4,IW)
      IF (.NOT.M%SOLID(IC)) M%WALL_INDEX_BACK(IW) = M%WALL_INDEX(IC,IOR)
      IF (     M%SOLID(IC)) THEN
         SELECT CASE(IOR)
            CASE(-1)
               II=II+1
            CASE( 1)
               II=II-1
            CASE(-2)
               JJ=JJ+1
            CASE( 2)
               JJ=JJ-1
            CASE(-3)
               KK=KK+1
            CASE( 3)
               KK=KK-1
         END SELECT
         IC = M%CELL_INDEX(II,JJ,KK)
         M%WALL_INDEX_BACK(IW) = M%WALL_INDEX(IC,IOR)
      ENDIF
   ENDIF
ENDDO

! Set clocks and counters related to frequency of solid phase conduction updates

M%BC_CLOCK     = T_BEGIN
M%WALL_COUNTER = 0

! Set clock for boudary fuel vegetation model

M%VEG_CLOCK_BC = T_BEGIN
 
! Allocate arrays for storing velocity boundary condition info
 
N_EDGES_DIM = 4*(IBP1*JBP1+IBP1*KBP1+JBP1*KBP1)
IF (EVACUATION_ONLY(NM)) N_EDGES_DIM = 4*(IBP1*KBP1+JBP1*KBP1)
DO N=1,M%N_OBST
   OB=>M%OBSTRUCTION(N)
   IPTS = OB%I2-OB%I1
   JPTS = OB%J2-OB%J1
   KPTS = OB%K2-OB%K1
   IF (EVACUATION_ONLY(NM)) THEN
      N_EDGES_DIM = N_EDGES_DIM + 4*(IPTS*KPTS+JPTS*KPTS)
   ELSE
      N_EDGES_DIM = N_EDGES_DIM + 4*(IPTS*JPTS+IPTS*KPTS+JPTS*KPTS)
   ENDIF
ENDDO

ALLOCATE(M%IJKE(16,N_EDGES_DIM),STAT=IZERO)
CALL ChkMemErr('INIT','IJKE',IZERO)   
M%IJKE  = 0
ALLOCATE(M%OME_E(0:N_EDGES_DIM,-2:2),STAT=IZERO)
CALL ChkMemErr('INIT','OME_E',IZERO)  
M%OME_E = 0._EB
ALLOCATE(M%TAU_E(0:N_EDGES_DIM,-2:2),STAT=IZERO)
CALL ChkMemErr('INIT','TAU_E',IZERO)  
M%TAU_E = 0._EB
ALLOCATE(M%EDGE_TYPE(N_EDGES_DIM,2),STAT=IZERO)
CALL ChkMemErr('INIT','EDGE_TYPE',IZERO)  
M%EDGE_TYPE = SOLID_EDGE
ALLOCATE(M%EDGE_INTERPOLATION_FACTOR(N_EDGES_DIM,2),STAT=IZERO)
CALL ChkMemErr('INIT','EDGE_INTERPOLATION_FACTOR',IZERO)  
M%EDGE_INTERPOLATION_FACTOR = 1._EB
 
! Initialize and allocate lagrangian particle/droplet arrays

M%NLP = 0
M%NLPDIM = 1000
IF (DROPLET_FILE .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   ALLOCATE(M%DROPLET(M%NLPDIM),STAT=IZERO)
   CALL ChkMemErr('INIT','DROPLET',IZERO)
ENDIF
 
! Allocate array to hold character strings for Smokeview file
 
M%N_STRINGS     =   0
M%N_STRINGS_MAX = 100
ALLOCATE(M%STRING(M%N_STRINGS_MAX),STAT=IZERO)
CALL ChkMemErr('INIT','STRING',IZERO)
 
! Set up arrays to hold velocity boundary condition info
 
CALL INITIALIZE_EDGES
 
! Initialize Pressure solver
 
CALL INITIALIZE_POISSON_SOLVER
IF (IERR/=0) RETURN
 
! Determine which wall cells to assign for solid phase thermocouples and profiles
 
CALL INITIALIZE_DEVC
IF (IERR/=0) RETURN

CALL INITIALIZE_PROF
IF (IERR/=0) RETURN
 
! Initialize Mesh Exchange
 
CALL INITIALIZE_INTERPOLATION
 
CONTAINS
 
 
SUBROUTINE INITIALIZE_EDGES
 
! Set up edge arrays for velocity boundary conditions
 
INTEGER I,J,K,N
 
CALL POINT_TO_MESH(NM)
 
N_EDGES = 0
 
! Arguments for DEFINE_EDGE(I,J,K,IOR,IEC,NM,I_OBST)
 
DO K=0,KBAR
   DO J=0,JBAR
      IF (J>0) CALL DEFINE_EDGE(   0,J,K, 1,2,NM,0,IERR)
      IF (J>0) CALL DEFINE_EDGE(IBAR,J,K,-1,2,NM,0,IERR)
      IF (K>0) CALL DEFINE_EDGE(   0,J,K, 1,3,NM,0,IERR)
      IF (K>0) CALL DEFINE_EDGE(IBAR,J,K,-1,3,NM,0,IERR)
   ENDDO
ENDDO
DO K=0,KBAR
   DO I=0,IBAR
      IF (I>0) CALL DEFINE_EDGE(I,   0,K, 2,1,NM,0,IERR)
      IF (I>0) CALL DEFINE_EDGE(I,JBAR,K,-2,1,NM,0,IERR)
      IF (K>0) CALL DEFINE_EDGE(I,   0,K, 2,3,NM,0,IERR)
      IF (K>0) CALL DEFINE_EDGE(I,JBAR,K,-2,3,NM,0,IERR)
   ENDDO
ENDDO
DO J=0,JBAR
   DO I=0,IBAR
      IF (I>0) CALL DEFINE_EDGE(I,J,   0, 3,1,NM,0,IERR)
      IF (I>0) CALL DEFINE_EDGE(I,J,KBAR,-3,1,NM,0,IERR)
      IF (J>0) CALL DEFINE_EDGE(I,J,   0, 3,2,NM,0,IERR)
      IF (J>0) CALL DEFINE_EDGE(I,J,KBAR,-3,2,NM,0,IERR)
   ENDDO
ENDDO

IF (IERR/=0) RETURN

OBST_LOOP_3: DO N=1,N_OBST
   OB => OBSTRUCTION(N)
   DO K=OB%K1,OB%K2
      DO J=OB%J1,OB%J2
         IF (J>OB%J1) CALL DEFINE_EDGE(OB%I1,J,K,-1,2,NM,N,IERR)
         IF (J>OB%J1) CALL DEFINE_EDGE(OB%I2,J,K, 1,2,NM,N,IERR)
         IF (K>OB%K1) CALL DEFINE_EDGE(OB%I1,J,K,-1,3,NM,N,IERR)
         IF (K>OB%K1) CALL DEFINE_EDGE(OB%I2,J,K, 1,3,NM,N,IERR)
      ENDDO
   ENDDO
   DO K=OB%K1,OB%K2
      DO I=OB%I1,OB%I2
         IF (I>OB%I1) CALL DEFINE_EDGE(I,OB%J1,K,-2,1,NM,N,IERR)
         IF (I>OB%I1) CALL DEFINE_EDGE(I,OB%J2,K, 2,1,NM,N,IERR)
         IF (K>OB%K1) CALL DEFINE_EDGE(I,OB%J1,K,-2,3,NM,N,IERR)
         IF (K>OB%K1) CALL DEFINE_EDGE(I,OB%J2,K, 2,3,NM,N,IERR)
      ENDDO
   ENDDO
   DO J=OB%J1,OB%J2
      DO I=OB%I1,OB%I2
         IF (I>OB%I1) CALL DEFINE_EDGE(I,J,OB%K1,-3,1,NM,N,IERR)
         IF (I>OB%I1) CALL DEFINE_EDGE(I,J,OB%K2, 3,1,NM,N,IERR)
         IF (J>OB%J1) CALL DEFINE_EDGE(I,J,OB%K1,-3,2,NM,N,IERR)
         IF (J>OB%J1) CALL DEFINE_EDGE(I,J,OB%K2, 3,2,NM,N,IERR)
      ENDDO
   ENDDO
ENDDO OBST_LOOP_3

IF (N_EDGES>N_EDGES_DIM .AND. EVACUATION_ONLY(NM)) THEN
   WRITE(LU_ERR,'(A,I2,A,2I8)') 'ERROR: Edges memory; Mesh: ',NM,', n_edges, n_edges_dim ',N_EDGES, N_EDGES_DIM
   PROCESS_STOP_STATUS = SETUP_STOP
   IERR = 1
ENDIF
 
END SUBROUTINE INITIALIZE_EDGES
 
 
SUBROUTINE INITIALIZE_POISSON_SOLVER

USE POIS, ONLY: H3CZIS,H2CZIS,H3CSIS,H2CYIS
REAL(EB) :: XLM,XMU
INTEGER  :: N
INTEGER, POINTER :: ITRN,JTRN,KTRN,LBC,MBC,NBC
INTEGER, POINTER, DIMENSION(:) :: NOC
TYPE (VENTS_TYPE), POINTER :: VT
 
! Allocate major arrays
 
ITRN =>M%ITRN 
JTRN =>M%JTRN
KTRN =>M%KTRN
LBC =>M%LBC
MBC =>M%MBC
NBC =>M%NBC
NOC=>TRANS(NM)%NOC
IF (NOC(1)==0 .AND. NOC(2)==0 .AND. NOC(3)==0) M%IPS=0
IF (NOC(1)/=0 .AND. NOC(2)==0 .AND. NOC(3)==0) M%IPS=1
IF (NOC(1)==0 .AND. NOC(2)/=0 .AND. NOC(3)==0) M%IPS=2
IF (NOC(1)==0 .AND. NOC(2)==0 .AND. NOC(3)/=0) M%IPS=3
IF (NOC(1)/=0 .AND. NOC(2)/=0 .AND. NOC(3)==0) M%IPS=4
IF (NOC(1)/=0 .AND. NOC(2)==0 .AND. NOC(3)/=0) M%IPS=5
IF (NOC(1)==0 .AND. NOC(2)/=0 .AND. NOC(3)/=0) M%IPS=6
IF (EVACUATION_ONLY(NM)                      ) M%IPS=7
IF (NOC(1)/=0 .AND. NOC(2)/=0 .AND. NOC(3)/=0) THEN
   WRITE(LU_ERR,'(A,I3)') 'ERROR: Stretch at most 2 coordinate directions in MESH ',NM
   PROCESS_STOP_STATUS = SETUP_STOP
   IERR = 1
   RETURN
ENDIF
 
IF (M%IPS<=1 .OR. M%IPS==4) THEN
   ITRN = IBP1
   IF (JBAR>1) JTRN = JBP1
   IF (JBAR==1) JTRN = 1
   KTRN = KBP1
   
   ! pressure periodic boundary conditions
   IF (FISHPAK_BC(1)==0) ITRN=IBAR
   IF (FISHPAK_BC(2)==0) JTRN=JBAR
   IF (FISHPAK_BC(3)==0) KTRN=KBAR
ENDIF
 
IF (M%IPS==2) THEN
   ITRN = JBP1
   JTRN = IBP1
   KTRN = KBP1
   ALLOCATE(M%BZST(JBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZST',IZERO)
   ALLOCATE(M%BZFT(JBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZFT',IZERO)
ENDIF
 
IF (M%IPS==3 .OR. M%IPS==6) THEN
   ITRN = KBP1
   IF (JBAR>1) JTRN = JBP1
   IF (JBAR==1) JTRN = 1
   KTRN = IBP1
   ALLOCATE(M%BXST(KBP1,JTRN),STAT=IZERO)
   CALL ChkMemErr('INIT','BXST',IZERO)
   ALLOCATE(M%BXFT(KBP1,JTRN),STAT=IZERO)
   CALL ChkMemErr('INIT','BXFT',IZERO)
   ALLOCATE(M%BYST(KBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BYST',IZERO)
   ALLOCATE(M%BYFT(KBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BYFT',IZERO)
   ALLOCATE(M%BZST(JTRN,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZST',IZERO)
   ALLOCATE(M%BZFT(JTRN,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZFT',IZERO)
ENDIF
 
IF (M%IPS==5) THEN
   ITRN = IBP1
   JTRN = KBP1
   KTRN = JBP1
   ALLOCATE(M%BXST(KBP1,JBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BXST',IZERO)
   ALLOCATE(M%BXFT(KBP1,JBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BXFT',IZERO)
ENDIF
 
IF (M%IPS==7) THEN
   ITRN = IBP1
   JTRN = JBP1
   KTRN = 1
ENDIF
 
IF (M%IPS<=3 .OR. M%IPS==7) THEN
   M%LSAVE = (ITRN+1)*JTRN*KTRN+7*ITRN+5*JTRN+6*KTRN+56
   M%LWORK = (ITRN+1)*JTRN*KTRN
ELSE
   N_LOOP: DO N=1,50
      IF ((JTRN+1)<=2**N) EXIT N_LOOP
   ENDDO N_LOOP
   M%LSAVE = KTRN*(6*N*(2**N)+2*N+19)+8*ITRN+7*JTRN+38
   M%LWORK = JTRN*(ITRN*(KTRN+1)+1)
ENDIF
 
ALLOCATE(M%SAVE1(-3:M%LSAVE),STAT=IZERO)
CALL ChkMemErr('INIT','SAVE1',IZERO)
IF (EMBEDDED_MESH) THEN
   ALLOCATE(M%SAVE2(-3:M%LSAVE),STAT=IZERO)
   CALL ChkMemErr('INIT','SAVE2',IZERO)
ENDIF
ALLOCATE(M%WORK(M%LWORK),STAT=IZERO)
CALL ChkMemErr('INIT','WORK',IZERO)
ALLOCATE(M%PRHS(ITRN,JTRN,KTRN),STAT=IZERO)
CALL ChkMemErr('INIT','PRHS',IZERO)
IF (KBAR>1) THEN
   IF (JBAR>1) ALLOCATE(M%BXS(JBP1,KBP1),STAT=IZERO)
   IF (JBAR==1) ALLOCATE(M%BXS(1,KBP1)   ,STAT=IZERO)
ELSE
               ALLOCATE(M%BXS(JBP1,1)   ,STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BXS',IZERO)
IF (KBAR>1) THEN
   IF (JBAR>1) ALLOCATE(M%BXF(JBP1,KBP1),STAT=IZERO)
   IF (JBAR==1) ALLOCATE(M%BXF(1,KBP1)   ,STAT=IZERO)
ELSE
               ALLOCATE(M%BXF(JBP1,1)   ,STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BXF',IZERO)
IF (KBAR>1) THEN
   ALLOCATE(M%BYS(IBP1,KBP1),STAT=IZERO)
ELSE
   ALLOCATE(M%BYS(IBP1,1),STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BYS',IZERO)
IF (KBAR>1) THEN
   ALLOCATE(M%BYF(IBP1,KBP1),STAT=IZERO)
ELSE
   ALLOCATE(M%BYF(IBP1,1),STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BYF',IZERO)
IF (JBAR>1) ALLOCATE(M%BZS(IBP1,JBP1),STAT=IZERO)
IF (JBAR==1) ALLOCATE(M%BZS(IBP1,1)   ,STAT=IZERO)
CALL ChkMemErr('INIT','BZS',IZERO)
IF (JBAR>1) ALLOCATE(M%BZF(IBP1,JBP1),STAT=IZERO)
IF (JBAR==1) ALLOCATE(M%BZF(IBP1,1)   ,STAT=IZERO)
CALL ChkMemErr('INIT','BZF',IZERO)
 
M%SAVE1  = 0._EB
IF (EMBEDDED_MESH) M%SAVE2 = 0._EB
M%WORK  = 0._EB
M%PRHS  = 0._EB
M%BXS   = 0._EB
M%BXF   = 0._EB
M%BYS   = 0._EB
M%BYF   = 0._EB
M%BZS   = 0._EB
M%BZF   = 0._EB
 
! Initialize pressure solver   
 
XLM = 0._EB         ! No Helmholtz equation
XMU = 0._EB         ! No Helmholtz equation
LBC = 3
MBC = 3
NBC = 3

! Look for OPEN vents -- this will change the entire face to DIRICHLET BCs
 
VENT_LOOP: DO N=1,M%N_VENT
   VT => M%VENTS(N)
   IF (VT%BOUNDARY_TYPE /= OPEN_BOUNDARY) CYCLE VENT_LOOP
   IF (VT%I1==0 .AND. VT%I2==0) THEN
      IF (LBC==3) LBC = 2
      IF (LBC==4) LBC = 1
   ENDIF
   IF (VT%I1==M%IBAR .AND. VT%I2==M%IBAR) THEN
      IF (LBC==3) LBC = 4
      IF (LBC==2) LBC = 1
   ENDIF
   IF (VT%J1==0 .AND. VT%J2==0) THEN
      IF (MBC==3) MBC = 2
      IF (MBC==4) MBC = 1
   ENDIF
   IF (VT%J1==M%JBAR .AND. VT%J2==M%JBAR) THEN
      IF (MBC==3) MBC = 4
      IF (MBC==2) MBC = 1
   ENDIF
   IF (VT%K1==0 .AND. VT%K2==0) THEN
      IF (NBC==3) NBC = 2
      IF (NBC==4) NBC = 1
   ENDIF
   IF (VT%K1==M%KBAR .AND. VT%K2==M%KBAR) THEN
      IF (NBC==3) NBC = 4
      IF (NBC==2) NBC = 1
   ENDIF
ENDDO VENT_LOOP
 
DO IW=1,M%NEWC
   IF (M%IJKW(9,IW)>0) THEN
      SELECT CASE(M%IJKW(4,IW))
         CASE( 1)
            IF (LBC==3) LBC = 2
            IF (LBC==4) LBC = 1
         CASE(-1)
            IF (LBC==3) LBC = 4
            IF (LBC==2) LBC = 1
         CASE( 2)
            IF (MBC==3) MBC = 2
            IF (MBC==4) MBC = 1
         CASE(-2)
            IF (MBC==3) MBC = 4
            IF (MBC==2) MBC = 1
         CASE( 3)
            IF (NBC==3) NBC = 2
            IF (NBC==4) NBC = 1
         CASE(-3)
            IF (NBC==3) NBC = 4
            IF (NBC==2) NBC = 1
      END SELECT
   ENDIF
ENDDO

! Periodic pressure boundary conditions for CrayFishpak

IF (FISHPAK_BC(1)==0) LBC=0
IF (FISHPAK_BC(2)==0) MBC=0
IF (FISHPAK_BC(3)==0) NBC=0
 
! Poisson solver with stretching in the 1st coordinate
 
SELECT_POISSON_SOLVER: SELECT CASE(M%IPS)

   CASE (0:1) SELECT_POISSON_SOLVER
      IF (.NOT.TWO_D) CALL H3CZIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,ZS,ZF,KBAR,NBC,M%HX,XLM,ITRN,JTRN,IERR,M%SAVE1)
      IF (TWO_D .AND. .NOT.CYLINDRICAL) CALL H2CZIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,M%HX,XLM,ITRN,IERR,M%SAVE1)
      IF (TWO_D .AND. CYLINDRICAL) THEN
         IF (XS==0._EB .AND. LBC==1) LBC = 5
         IF (XS==0._EB .AND. LBC==2) LBC = 6
         IF (XS==0._EB .AND. LBC==3) LBC = 6
         IF (XS==0._EB .AND. LBC==4) LBC = 5
         CALL H2CYIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,XLM,XMU,ITRN,IERR,M%SAVE1)
      ENDIF
      IF (EMBEDDED_MESH) THEN
         IF (.NOT.TWO_D) CALL H3CZIS(XS,XF,IBAR,1,YS,YF,JBAR,1,ZS,ZF,KBAR,1,M%HX,XLM,ITRN,JTRN,IERR,M%SAVE2)
         IF (TWO_D .AND. .NOT.CYLINDRICAL) CALL H2CZIS(XS,XF,IBAR,1,ZS,ZF,KBAR,1,M%HX,XLM,ITRN,IERR,M%SAVE2)
      ENDIF
   CASE (2) SELECT_POISSON_SOLVER
      CALL H3CZIS(YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,M%HY,XLM,ITRN,JTRN,IERR,M%SAVE1)
   CASE (3) SELECT_POISSON_SOLVER
      IF (TWO_D) THEN
         CALL H2CZIS(ZS,ZF,KBAR,NBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,IERR,M%SAVE1)
      ELSE
         CALL H3CZIS(ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,JTRN,IERR,M%SAVE1)
      ENDIF
   CASE (4) SELECT_POISSON_SOLVER
      CALL H3CSIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,ZS,ZF,KBAR,NBC,XLM,ITRN,JTRN,IERR,M%SAVE1,M%WORK,M%HX,M%HY)
   CASE (5) SELECT_POISSON_SOLVER
      IF (TWO_D) THEN
         CALL H2CZIS(ZS,ZF,KBAR,NBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,IERR,M%SAVE1)
      ELSE
         CALL H3CSIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XLM,ITRN,JTRN,IERR,M%SAVE1,M%WORK,M%HX,M%HZ)
      ENDIF
   CASE (6) SELECT_POISSON_SOLVER
      CALL H3CSIS(ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,XLM,ITRN,JTRN,IERR,M%SAVE1,M%WORK,M%HZ,M%HY)
   CASE (7) SELECT_POISSON_SOLVER
      CALL H2CZIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,M%HX,XLM,ITRN,IERR,M%SAVE1)

END SELECT SELECT_POISSON_SOLVER
 
! Specify the pressure boundary condition for each wall cell

WALL_CELL_LOOP: DO IW=1,M%NEWC
   IOR = M%IJKW(4,IW)
   SELECT CASE(IOR)
      CASE( 1)
         IF (LBC==3 .OR. LBC==4 .OR. LBC==6) M%PRESSURE_BC_INDEX(IW) = NEUMANN
         IF (LBC==1 .OR. LBC==2 .OR. LBC==5) M%PRESSURE_BC_INDEX(IW) = DIRICHLET
      CASE(-1)
         IF (LBC==2 .OR. LBC==3 .OR. LBC==6) M%PRESSURE_BC_INDEX(IW) = NEUMANN
         IF (LBC==1 .OR. LBC==4 .OR. LBC==5) M%PRESSURE_BC_INDEX(IW) = DIRICHLET
      CASE( 2)
         IF (MBC==3 .OR. MBC==4) M%PRESSURE_BC_INDEX(IW) = NEUMANN
         IF (MBC==1 .OR. MBC==2) M%PRESSURE_BC_INDEX(IW) = DIRICHLET
      CASE(-2)
         IF (MBC==3 .OR. MBC==2) M%PRESSURE_BC_INDEX(IW) = NEUMANN
         IF (MBC==1 .OR. MBC==4) M%PRESSURE_BC_INDEX(IW) = DIRICHLET
      CASE( 3)
         IF (NBC==3 .OR. NBC==4) M%PRESSURE_BC_INDEX(IW) = NEUMANN
         IF (NBC==1 .OR. NBC==2) M%PRESSURE_BC_INDEX(IW) = DIRICHLET
      CASE(-3)
         IF (NBC==3 .OR. NBC==2) M%PRESSURE_BC_INDEX(IW) = NEUMANN
         IF (NBC==1 .OR. NBC==4) M%PRESSURE_BC_INDEX(IW) = DIRICHLET
   END SELECT
ENDDO WALL_CELL_LOOP

! Check for errors with Poisson solver initialization
 
IF (IERR/=0) THEN
   WRITE(LU_ERR,'(A,I2,A,I3)') 'ERROR: Poisson initialization error, Number=',IERR, ', Mesh=',NM
   PROCESS_STOP_STATUS = SETUP_STOP
   RETURN
ENDIF
 
END SUBROUTINE INITIALIZE_POISSON_SOLVER
 
 
 
SUBROUTINE INITIALIZE_DEVC

INTEGER :: III
TYPE (DEVICE_TYPE), POINTER :: DV
 
DEVICE_LOOP: DO N=1,N_DEVC
   DV => DEVICE(N)
   IF (NM/=DV%MESH) CYCLE DEVICE_LOOP

   SOLID_OR_GAS: IF (DV%OUTPUT_INDEX<0)  THEN

      II  = GINV(DV%X-M%XS,1,NM)*M%RDXI   + 1._EB
      JJ  = GINV(DV%Y-M%YS,2,NM)*M%RDETA  + 1._EB
      KK  = GINV(DV%Z-M%ZS,3,NM)*M%RDZETA + 1._EB
      IOR = DV%IOR
      CALL GET_IW(II,JJ,KK,IOR,IW)
      IF (IW>0 .OR. DV%STATISTICS/='null') THEN
         DV%IW = IW
         IF (OUTPUT_QUANTITY(DV%OUTPUT_INDEX)%INSIDE_SOLID) THEN
            IF (IW>0) THEN
               IBC = M%IJKW(5,IW)
            ELSE
               IBC = DV%SURF_INDEX
            ENDIF
            IF (SURFACE(IBC)%THERMAL_BC_INDEX /= THERMALLY_THICK) THEN
               WRITE(LU_ERR,'(A,I3,A)') 'ERROR: DEViCe ',N, ' must be associated with a heat-conducting surface'
               PROCESS_STOP_STATUS = SETUP_STOP
               IERR = 1
               RETURN
            ENDIF
            DV%I_DEPTH = SURFACE(IBC)%N_CELLS
            DO III=SURFACE(IBC)%N_CELLS,1,-1
               IF (DV%DEPTH<=SURFACE(IBC)%X_S(III)) DV%I_DEPTH = III
            ENDDO
         ENDIF
      ELSE
         WRITE(LU_ERR,'(A,I4,A)') 'ERROR: Reposition DEVC No.',DV%ORDINAL, '. FDS cannot determine which boundary cell to assign'
         PROCESS_STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF

   ELSE SOLID_OR_GAS

      IF (DV%VIRTUAL_WALL_INDEX<1) CYCLE DEVICE_LOOP
      IW  = DV%VIRTUAL_WALL_INDEX
      IBC = M%IJKW(5,IW)
      DV%I_DEPTH = SURFACE(IBC)%N_CELLS
      DO III=SURFACE(IBC)%N_CELLS,1,-1
         IF (DV%DEPTH<=SURFACE(IBC)%X_S(III)) DV%I_DEPTH = III
      ENDDO

   ENDIF SOLID_OR_GAS

ENDDO DEVICE_LOOP
 
END SUBROUTINE INITIALIZE_DEVC
 
 
SUBROUTINE INITIALIZE_PROF

INTEGER :: NN
LOGICAL :: SUCCESS
TYPE (PROFILE_TYPE), POINTER :: PF

IF (EVACUATION_ONLY(NM)) RETURN
PROF_LOOP: DO N=1,N_PROF
   PF => PROFILE(N)
   IF (NM/=PF%MESH) CYCLE PROF_LOOP
   II  = GINV(PF%X-M%XS,1,NM)*M%RDXI   + 1._EB
   JJ  = GINV(PF%Y-M%YS,2,NM)*M%RDETA  + 1._EB
   KK  = GINV(PF%Z-M%ZS,3,NM)*M%RDZETA + 1._EB
   IOR = PF%IOR
   CALL GET_IW(II,JJ,KK,IOR,IW)
   IF (IW>0) THEN
      PF%IW = IW
      SF => SURFACE(IJKW(5,IW))
      IF (.NOT.SF%THERMALLY_THICK) THEN
         WRITE(LU_ERR,'(A,I3,A)') 'ERROR: PROFile ',N, ' must be associated with a heat-conducting surface'
         PROCESS_STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF
      IF (PF%QUANTITY /= 'TEMPERATURE' .AND. PF%QUANTITY /= 'DENSITY') THEN
         SUCCESS = .FALSE.
         DO NN=1,SF%N_MATL
            IF (PF%QUANTITY==SF%MATL_NAME(NN)) SUCCESS = .TRUE.
         ENDDO
         IF (.NOT.SUCCESS) THEN
            WRITE(LU_ERR,'(A,A,A)') 'ERROR: QUANTITY ',TRIM(PF%QUANTITY), ' is not appropriate for the designated location'
            PROCESS_STOP_STATUS = SETUP_STOP
            IERR = 1
            RETURN
         ENDIF
      ENDIF
   ELSE
      WRITE(LU_ERR,'(A,I4,A)') 'ERROR: Reposition PROF No.',PF%ORDINAL, '. FDS cannot determine which boundary cell to assign'
      PROCESS_STOP_STATUS = SETUP_STOP
      IERR = 1
      RETURN
   ENDIF
ENDDO PROF_LOOP
 
END SUBROUTINE INITIALIZE_PROF


SUBROUTINE GET_IW(II,JJ,KK,IOR,IW)

INTEGER :: II,JJ,KK,IOR,IW

IC  = M%CELL_INDEX(II,JJ,KK)
 
IF (M%SOLID(IC)) THEN
   SELECT CASE(IOR)
      CASE(-1)
         IF (II>0)      II = II-1
      CASE( 1)
         IF (II<M%IBP1) II = II+1
      CASE(-2)
         IF (JJ>0)      JJ = JJ-1
      CASE( 2)
         IF (JJ<M%JBP1) JJ = JJ+1
      CASE(-3)
         IF (KK>0)      KK = KK-1
      CASE( 3)
         IF (KK<M%KBP1) KK = KK+1
   END SELECT
ENDIF
 
IC  = M%CELL_INDEX(II,JJ,KK)
IW  = M%WALL_INDEX(IC,-IOR)
 
IF (IW<=0) THEN
   SELECT CASE(IOR)
      CASE(-1)
         IF (II>0)      IC = M%CELL_INDEX(II-1,JJ,KK)
      CASE( 1)
         IF (II<M%IBP1) IC = M%CELL_INDEX(II+1,JJ,KK)
      CASE(-2)
         IF (JJ>0)      IC = M%CELL_INDEX(II,JJ-1,KK)
      CASE( 2)
         IF (JJ<M%JBP1) IC = M%CELL_INDEX(II,JJ+1,KK)
      CASE(-3)
         IF (KK>0)      IC = M%CELL_INDEX(II,JJ,KK-1)
      CASE( 3)
         IF (KK<M%KBP1) IC = M%CELL_INDEX(II,JJ,KK+1)
   END SELECT
   IW = M%WALL_INDEX(IC,-IOR)
ENDIF

END SUBROUTINE GET_IW



SUBROUTINE INITIALIZE_INTERPOLATION
 
! Create arrays by which info is to exchanged across meshes
 
INTEGER :: NOM,I,J,K
TYPE (MESH_TYPE), POINTER :: M2
 
IF (NM==1) RETURN
IF (EVACUATION_ONLY(NM)) RETURN
 
ALLOCATE(M%INTERPOLATED_MESH(1:M%IBAR,1:M%JBAR,1:M%KBAR),  STAT=IZERO)
CALL ChkMemErr('INIT','INTERPOLATED_MESH',IZERO)  
M%INTERPOLATED_MESH = 0
 
DO K=1,M%KBAR
   DO J=1,M%JBAR
      DO I=1,M%IBAR
         OTHER_MESH_LOOP: DO NOM=1,NM-1
            M2=>MESHES(NOM)
            IF (M%X(I-1)>=M2%XS .AND. M%X(I)<=M2%XF .AND.  M%Y(J-1)>=M2%YS .AND. M%Y(J)<=M2%YF .AND. &
                M%Z(K-1)>=M2%ZS .AND. M%Z(K)<=M2%ZF) THEN 
               M%INTERPOLATED_MESH(I,J,K) = NOM
               EXIT OTHER_MESH_LOOP
            ENDIF
         ENDDO OTHER_MESH_LOOP
      ENDDO
   ENDDO
ENDDO
 
END SUBROUTINE INITIALIZE_INTERPOLATION
 
END SUBROUTINE INITIALIZE_MESH_VARIABLES
 
 
 
SUBROUTINE INITIALIZE_GLOBAL_VARIABLES
USE CONTROL_VARIABLES, ONLY: N_CTRL

! Initialize time, printout and plot clocks
 
ALLOCATE(PART_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PART_CLOCK',IZERO) 
ALLOCATE(ISOF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','ISOF_CLOCK',IZERO) 
ALLOCATE(BNDF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','BNDF_CLOCK',IZERO) 
ALLOCATE(SLCF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','SLCF_CLOCK',IZERO) 
ALLOCATE(CORE_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','CORE_CLOCK',IZERO) 
ALLOCATE(PL3D_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PL3D_CLOCK',IZERO) 
ALLOCATE(PROF_CLOCK(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PROF_CLOCK',IZERO) 
 
ICYC              = 0
PART_CLOCK        = T_BEGIN
DEVC_CLOCK        = T_BEGIN      
CTRL_CLOCK        = T_BEGIN      
PROF_CLOCK        = T_BEGIN      
PL3D_CLOCK        = T_BEGIN + DT_PL3D
ISOF_CLOCK        = T_BEGIN      
SLCF_CLOCK        = T_BEGIN   
BNDF_CLOCK        = T_BEGIN     
CORE_CLOCK        = T_BEGIN + DT_RESTART
FLUSH_CLOCK       = T_BEGIN + DT_FLUSH
HRR_CLOCK         = T_BEGIN
MINT_CLOCK        = T_BEGIN      
VEG_CLOCK         = T_BEGIN
IF (PERIODIC_TEST==2) SPEC_CLOCK = T_BEGIN

IF (N_DEVC==0)      DEVC_CLOCK = 1.E10_EB
IF (N_CTRL==0)      CTRL_CLOCK = 1.E10_EB
IF (N_PROF==0)      PROF_CLOCK = 1.E10_EB
IF (N_ISOF==0)      ISOF_CLOCK = 1.E10_EB
IF (N_BNDF==0)      BNDF_CLOCK = 1.E10_EB
IF (.NOT.MASS_FILE) MINT_CLOCK = 1.E10_EB

ALLOCATE(HRR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','HRR',IZERO) 
HRR = 0._EB
ALLOCATE(RHRR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','RHRR',IZERO) 
RHRR = 0._EB
ALLOCATE(CHRR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','CHRR',IZERO) 
CHRR = 0._EB
ALLOCATE(FHRR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','FHRR',IZERO) 
FHRR = 0._EB
ALLOCATE(MLR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MLR',IZERO) 
MLR = 0._EB
ALLOCATE(FUEL_DROPLET_MLR(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','FUEL_DROPLET_MLR',IZERO) 
FUEL_DROPLET_MLR = 0._EB
ALLOCATE(HRR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','HRR_SUM',IZERO)  
HRR_SUM=0._EB
ALLOCATE(RHRR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','RHRR_SUM',IZERO) 
RHRR_SUM=0._EB
ALLOCATE(CHRR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','CHRR_SUM',IZERO) 
CHRR_SUM=0._EB
ALLOCATE(FHRR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','FHRR_SUM',IZERO) 
FHRR_SUM=0._EB
ALLOCATE(MLR_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MLR_SUM',IZERO)  
MLR_SUM=0._EB
ALLOCATE(HRR_TIME_INTERVAL(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MLR_COUNT',IZERO)  
HRR_TIME_INTERVAL=0._EB
ALLOCATE(MINT(0:MAX_SPECIES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MINT',IZERO) 
ALLOCATE(MINT_SUM(0:MAX_SPECIES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MINT_SUM',IZERO) 
MINT_SUM=0._EB
ALLOCATE(MINT_TIME_INTERVAL(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MINT_TIME_INTERVAL',IZERO) 
MINT_TIME_INTERVAL=0._EB
ALLOCATE(I_MIN(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','I_MIN',IZERO) 
I_MIN = -10
ALLOCATE(I_MAX(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','I_MAX',IZERO) 
I_MAX = -10
ALLOCATE(J_MIN(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','J_MIN',IZERO) 
J_MIN = -10
ALLOCATE(J_MAX(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','J_MAX',IZERO) 
J_MAX = -10
ALLOCATE(K_MIN(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','K_MIN',IZERO) 
K_MIN = -10
ALLOCATE(K_MAX(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','K_MAX',IZERO) 
K_MAX = -10
ALLOCATE(NIC(NMESHES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','NIC',IZERO) 
NIC = 0
ALLOCATE(T_PER_STEP(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','T_PER_STEP',IZERO) 
T_PER_STEP = 0._EB
ALLOCATE(T_ACCUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','T_ACCUM',IZERO) 
T_ACCUM = 0._EB
ALLOCATE(NTCYC(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','NTCYC',IZERO) 
NTCYC = 0
ALLOCATE(NCYC(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','NCYC',IZERO) 
NCYC = 0

ALLOCATE(DSUM(N_ZONE,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','DSUM',IZERO) 
DSUM = 0._EB
ALLOCATE(PSUM(N_ZONE,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PSUM',IZERO) 
PSUM = 0._EB
ALLOCATE(USUM(N_ZONE,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','USUM',IZERO) 
USUM = 0._EB
ALLOCATE(ASUM(N_ZONE,0:N_ZONE,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','ASUM',IZERO) 
ASUM = 0._EB
ALLOCATE(ACTUAL_LEAK_AREA(0:N_ZONE,0:N_ZONE),STAT=IZERO)
CALL ChkMemErr('INIT','ACTUAL_LEAK_AREA',IZERO) 
ACTUAL_LEAK_AREA = 0._EB
ALLOCATE(FDS_LEAK_AREA(0:N_ZONE,0:N_ZONE,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','FDS_LEAK_AREA',IZERO) 
FDS_LEAK_AREA = 0._EB
ALLOCATE(FDS_LEAK_AREA_RATIO(0:N_ZONE,0:N_ZONE,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','FDS_LEAK_AREA_RATIO',IZERO) 
FDS_LEAK_AREA_RATIO = 1._EB
ALLOCATE(VELOCITY_ERROR_MAX(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VELOCITY_ERROR_MAX',IZERO) 
VELOCITY_ERROR_MAX = 0._EB
ALLOCATE(VELOCITY_ERROR_MAX_I(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VELOCITY_ERROR_MAX_I',IZERO) 
ALLOCATE(VELOCITY_ERROR_MAX_J(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VELOCITY_ERROR_MAX_J',IZERO) 
ALLOCATE(VELOCITY_ERROR_MAX_K(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VELOCITY_ERROR_MAX_K',IZERO) 

! Vegetation arrays

ALLOCATE(VEG_TOTAL_DRY_MASS(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VEG_TOTAL_DRY_MASS',IZERO)  
VEG_TOTAL_DRY_MASS=0._EB
ALLOCATE(VEG_TOTAL_MOIST_MASS(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VEG_TOTAL_MOIST_MASS',IZERO)  
VEG_TOTAL_MOIST_MASS=0._EB

! PERM (particle advection scheme)
IF (.FALSE.) THEN
   ! Initialize Moore-Penrose pseudoinverse
   ALLOCATE(MATRIX_BPLUS(12,8),STAT=IZERO)
   CALL ChkMemErr('INIT','MATRIX_BPLUS',IZERO)  
   MATRIX_BPLUS(1,:)  = (/-7,-2,-2,-1, 7, 2, 2, 1/)
   MATRIX_BPLUS(2,:)  = (/-2,-1,-7,-2, 2, 1, 7, 2/)
   MATRIX_BPLUS(3,:)  = (/-2,-7,-1,-2, 2, 7, 1, 2/)
   MATRIX_BPLUS(4,:)  = (/-1,-2,-2,-7, 1, 2, 2, 7/)
   MATRIX_BPLUS(5,:)  = (/-7,-2, 7, 2,-2,-1, 2, 1/)
   MATRIX_BPLUS(6,:)  = (/-2,-1, 2, 1,-7,-2, 7, 2/)
   MATRIX_BPLUS(7,:)  = (/-2,-7, 2, 7,-1,-2, 1, 2/)
   MATRIX_BPLUS(8,:)  = (/-1,-2, 1, 2,-2,-7, 2, 7/)
   MATRIX_BPLUS(9,:)  = (/-7, 7,-2, 2,-2, 2,-1, 1/)
   MATRIX_BPLUS(10,:) = (/-2, 2,-7, 7,-1, 1,-2, 2/)
   MATRIX_BPLUS(11,:) = (/-2, 2,-1, 1,-7, 7,-2, 2/)
   MATRIX_BPLUS(12,:) = (/-1, 1,-2, 2,-2, 2,-7, 7/)
   MATRIX_BPLUS = MATRIX_BPLUS/12._EB
ENDIF
 
END SUBROUTINE INITIALIZE_GLOBAL_VARIABLES
 
 

SUBROUTINE INIT_WALL_CELL(NM,I,J,K,I_OBST,IW,IOR,IBC,IERR)

! Initialize wall cell variables at external and obstruction boundaries
 
USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_WALL_ARRAYS
USE GEOMETRY_FUNCTIONS, ONLY : SEARCH_OTHER_MESHES
INTEGER  :: NM,NOM,IBC,IOR,ITER,IIO_MIN,IIO_MAX,JJO_MIN,JJO_MAX,KKO_MIN,KKO_MAX,I_VENT
INTEGER, INTENT(OUT) :: IERR
REAL(EB) :: PX,PY,PZ,X1,X2,Y1,Y2,Z1,Z2,T_ACTIVATE,XIN,YIN,ZIN,DIST
INTEGER  :: N,I_OBST,I,J,K,IBCX,IIG,JJG,KKG,IW,IIO,JJO,KKO,IC,ICG,NOM_CHECK(0:1)
INTEGER :: NSLICE
LOGICAL :: VENT_FOUND,ALIGNED
TYPE (MESH_TYPE), POINTER :: MM
TYPE (OBSTRUCTION_TYPE), POINTER :: OBX
TYPE (VENTS_TYPE), POINTER :: VT
 
IERR = 0
M=>MESHES(NM)
 
! Determine if a VENT covers the surface
 
I_VENT = 0
IBCX   = IBC
VENT_FOUND = .FALSE.

VENT_SEARCH_LOOP: DO N=1,M%N_VENT

   VT => M%VENTS(N)
   IF (I_OBST>0) THEN
      IF (VT%BOUNDARY_TYPE==OPEN_BOUNDARY)       CYCLE VENT_SEARCH_LOOP
      IF (.NOT.M%OBSTRUCTION(I_OBST)%ALLOW_VENT) CYCLE VENT_SEARCH_LOOP
   ENDIF
   IF (VT%IOR/=IOR) CYCLE VENT_SEARCH_LOOP

   IF (ABS(IOR)==1) THEN
      IF (IOR== 1 .AND. I/=VT%I1  ) CYCLE VENT_SEARCH_LOOP
      IF (IOR==-1 .AND. I/=VT%I1+1) CYCLE VENT_SEARCH_LOOP
      IF (J<VT%J1+1 .OR. J>VT%J2)   CYCLE VENT_SEARCH_LOOP
      IF (K<VT%K1+1 .OR. K>VT%K2)   CYCLE VENT_SEARCH_LOOP
   ENDIF
   IF (ABS(IOR)==2) THEN
      IF (IOR== 2 .AND. J/=VT%J1  ) CYCLE VENT_SEARCH_LOOP
      IF (IOR==-2 .AND. J/=VT%J1+1) CYCLE VENT_SEARCH_LOOP
      IF (I<VT%I1+1 .OR. I>VT%I2)   CYCLE VENT_SEARCH_LOOP
      IF (K<VT%K1+1 .OR. K>VT%K2)   CYCLE VENT_SEARCH_LOOP
   ENDIF
   IF (ABS(IOR)==3) THEN
      IF (IOR== 3 .AND. K/=VT%K1  ) CYCLE VENT_SEARCH_LOOP
      IF (IOR==-3 .AND. K/=VT%K1+1) CYCLE VENT_SEARCH_LOOP
      IF (I<VT%I1+1 .OR. I>VT%I2)   CYCLE VENT_SEARCH_LOOP
      IF (J<VT%J1+1 .OR. J>VT%J2)   CYCLE VENT_SEARCH_LOOP
   ENDIF

   ! Check if there are over-lapping VENTs

   IF (VENT_FOUND) THEN
      WRITE(LU_ERR,'(A,I0,A,3(I0,1X),A,I0,A)') 'WARNING: Two VENTs overlap in MESH ',NM,', Cell ',I,J,K,'. VENT ', &
                                           VT%ORDINAL,' rejected for that cell'
      EXIT VENT_SEARCH_LOOP
   ENDIF

   VENT_FOUND = .TRUE.

   ! Reassign the SURF index to be that of the VENT

   I_VENT = N
   IBCX   = VT%IBC

ENDDO VENT_SEARCH_LOOP

! Compute boundary cell physical coords (XW,YW,ZW) and area (AW)

M%IJKW(1,IW) = I
M%IJKW(2,IW) = J
M%IJKW(3,IW) = K
M%IJKW(4,IW) = IOR
M%IJKW(5,IW) = IBCX
M%OBST_INDEX_W(IW) = I_OBST
 
IF (ABS(IOR)==1) THEN
   IF (IOR== 1) THEN
      M%XW(IW)     = M%X(I)
      M%IJKW(6,IW) = I+1
      M%RDN(IW)    = M%RDXN(I) 
      IF (.NOT.EVACUATION_ONLY(NM)) M%AW(IW)     = M%R(I)*M%DY(J)*M%DZ(K)
      M%UW(IW)     = -U0
   ENDIF
   IF (IOR==-1) THEN
      M%XW(IW)     = M%X(I-1)
      M%IJKW(6,IW) = I-1
      M%RDN(IW)    = M%RDXN(I-1) 
      IF (.NOT.EVACUATION_ONLY(NM)) M%AW(IW)     = M%R(I-1)*M%DY(J)*M%DZ(K)
      M%UW(IW)     = U0
   ENDIF
   M%IJKW(7,IW) = J
   M%IJKW(8,IW) = K
   M%YW(IW) = 0.5_EB*(M%Y(J)+M%Y(J-1))  
   M%ZW(IW) = 0.5_EB*(M%Z(K)+M%Z(K-1))
ENDIF
IF (ABS(IOR)==2) THEN
   IF (IOR== 2) THEN
      M%YW(IW)     = M%Y(J)
      M%IJKW(7,IW) = J+1
      M%RDN(IW)    = M%RDYN(J) 
      M%UW(IW)     = -V0
   ENDIF
   IF (IOR==-2) THEN
      M%YW(IW)     = M%Y(J-1)
      M%IJKW(7,IW) = J-1
      M%RDN(IW)    = M%RDYN(J-1) 
      M%UW(IW)     = V0
   ENDIF
   M%IJKW(6,IW) = I
   M%IJKW(8,IW) = K
   M%XW(IW) = 0.5_EB*(M%X(I)+M%X(I-1))
   M%ZW(IW) = 0.5_EB*(M%Z(K)+M%Z(K-1))
   IF (.NOT.EVACUATION_ONLY(NM)) M%AW(IW) = M%DX(I)*M%DZ(K)
ENDIF
IF (ABS(IOR)==3) THEN
   IF (IOR== 3) THEN
      M%ZW(IW)     = M%Z(K)
      M%IJKW(8,IW) = K+1
      M%RDN(IW)    = M%RDZN(K) 
      M%UW(IW)     = -W0
   ENDIF
   IF (IOR==-3) THEN
      M%ZW(IW)     = M%Z(K-1)
      M%IJKW(8,IW) = K-1
      M%RDN(IW)    = M%RDZN(K-1) 
      M%UW(IW)     = W0
   ENDIF
   M%IJKW(6,IW) = I
   M%IJKW(7,IW) = J
   M%XW(IW) = 0.5_EB*(M%X(I)+M%X(I-1))
   M%YW(IW) = 0.5_EB*(M%Y(J)+M%Y(J-1))
   IF (.NOT.EVACUATION_ONLY(NM)) M%AW(IW) = M%DX(I)*M%RC(I)*M%DY(J)
ENDIF
 
IF (.NOT.EVACUATION_ONLY(NM)) THEN 
IF (M%AW(IW)>0._EB) M%RAW(IW) = 1._EB/M%AW(IW)
ENDIF

IF (IOR==0) THEN
   M%IJKW(6,IW) = I
   M%IJKW(7,IW) = J
   M%IJKW(8,IW) = K
ENDIF

! Fill array containing K index of terrain following slice, used in dump slice

DO NSLICE = 1, M%N_TERRAIN_SLCF
   IF (IOR==3) M%K_AGL_SLICE(I,J,NSLICE) = M%K_AGL_SLICE(I,J,NSLICE) + M%IJKW(8,IW) 
ENDDO

! Fill arrays containing K index and physical height of terrain for use in level set firespread 
IF (VEG_LEVEL_SET .AND. IOR==3) THEN 
  M%LS_Z_TERRAIN(M%IJKW(6,IW),M%IJKW(7,IW)) = M%Z(M%IJKW(8,IW)-1)
ENDIF

! Do not assign normal velocities at boundaries of evacuation meshes

IF (EVACUATION_ONLY(NM)) M%UW(IW) = 0._EB

! Save the wall index

IIG = M%IJKW(6,IW)
JJG = M%IJKW(7,IW)
KKG = M%IJKW(8,IW)
IC  = M%CELL_INDEX(I  ,J  ,K  )
ICG = M%CELL_INDEX(IIG,JJG,KKG)
M%WALL_INDEX(ICG,-IOR) = IW

! Use BOUNDARY_TYPE to indicate whether the boundary cell is blocked or on an obstruction that is HIDDEN

IF (IW<=M%NEWC .AND. I_OBST==0) M%BOUNDARY_TYPE(IW) = SOLID_BOUNDARY

IF (I_OBST>0) THEN
   M%BOUNDARY_TYPE(IW) = SOLID_BOUNDARY
   IF (IW<=M%NEWC) M%SOLID(IC) = .TRUE.
ENDIF

IF (M%SOLID(ICG)) M%BOUNDARY_TYPE(IW) = NULL_BOUNDARY
IF (SURFACE(IBC)%POROUS .AND. M%BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) M%BOUNDARY_TYPE(IW) = POROUS_BOUNDARY
IF (IOR==0) M%BOUNDARY_TYPE(IW) = VIRTUAL_BOUNDARY

! Assign the ZONE number to all boundary cells

IF (.NOT.EVACUATION_ONLY(NM)) M%PRESSURE_ZONE_WALL(IW) = M%PRESSURE_ZONE(IIG,JJG,KKG)
 
! Check for neighboring meshes in a multiple mesh calculation
 
CHECK_MESHES: IF (IW<=M%NEWC .AND.  .NOT.EVACUATION_ONLY(NM)) THEN
 
   IIO_MIN =  1000000
   IIO_MAX = -1000000
   JJO_MIN =  1000000
   JJO_MAX = -1000000
   KKO_MIN =  1000000
   KKO_MAX = -1000000
   NOM_CHECK = 0
   DO ITER=0,1
      XIN = M%XW(IW)
      YIN = M%YW(IW)
      ZIN = M%ZW(IW)
      IF (IBCX==PERIODIC_SURF_INDEX .AND. IOR== 1) XIN = XF_MAX
      IF (IBCX==PERIODIC_SURF_INDEX .AND. IOR==-1) XIN = XS_MIN
      IF (IBCX==PERIODIC_SURF_INDEX .AND. IOR== 2) YIN = YF_MAX
      IF (IBCX==PERIODIC_SURF_INDEX .AND. IOR==-2) YIN = YS_MIN
      IF (IBCX==PERIODIC_SURF_INDEX .AND. IOR== 3) ZIN = ZF_MAX
      IF (IBCX==PERIODIC_SURF_INDEX .AND. IOR==-3) ZIN = ZS_MIN
      IF (ABS(IOR)/=1) XIN = M%XW(IW) + (ITER*0.95_EB-0.475_EB)*(M%X(I)-M%X(I-1))
      IF (ABS(IOR)/=2) YIN = M%YW(IW) + (ITER*0.95_EB-0.475_EB)*(M%Y(J)-M%Y(J-1))
      IF (ABS(IOR)/=3) ZIN = M%ZW(IW) + (ITER*0.95_EB-0.475_EB)*(M%Z(K)-M%Z(K-1))
      IF (IOR== 1) XIN = XIN - 0.05*M%DX(0)
      IF (IOR==-1) XIN = XIN + 0.05*M%DX(M%IBP1)
      IF (IOR== 2) YIN = YIN - 0.05*M%DY(0)
      IF (IOR==-2) YIN = YIN + 0.05*M%DY(M%JBP1)
      IF (IOR== 3) ZIN = ZIN - 0.05*M%DZ(0)
      IF (IOR==-3) ZIN = ZIN + 0.05*M%DZ(M%KBP1)
      CALL SEARCH_OTHER_MESHES(XIN,YIN,ZIN,NOM,IIO,JJO,KKO) 
      NOM_CHECK(ITER) = NOM
      IF (NOM/=0) THEN
         IIO_MIN = MIN(IIO_MIN,IIO)
         IIO_MAX = MAX(IIO_MAX,IIO)
         JJO_MIN = MIN(JJO_MIN,JJO)
         JJO_MAX = MAX(JJO_MAX,JJO)
         KKO_MIN = MIN(KKO_MIN,KKO)
         KKO_MAX = MAX(KKO_MAX,KKO)
      ENDIF
   ENDDO
  
   ! Check to see if the current interpolated cell face spans more than one other mesh

   IF (NOM_CHECK(0)/=NOM_CHECK(1)) THEN
      WRITE(LU_ERR,'(A,I3,A,I3)') 'ERROR: MESH ',NM,' is out of alignment with MESH ',MAXVAL(NOM_CHECK)
      PROCESS_STOP_STATUS = SETUP_STOP
      IERR = 1
      RETURN
   ENDIF

   FOUND_OTHER_MESH: IF (NOM>0) THEN 
      MM=>MESHES(NOM)
      ALIGNED = .TRUE.
      IF ( (ABS(IOR)==2 .OR. ABS(IOR)==3) .AND. MM%DX(IIO_MIN)<=M%DX(I) .AND. &
            ABS( ((IIO_MAX-IIO_MIN+1)*MM%DX(IIO_MIN)-M%DX(I)) / MM%DX(IIO_MIN))>0.05 ) ALIGNED = .FALSE.
      IF ( (ABS(IOR)==1 .OR. ABS(IOR)==3) .AND. MM%DY(JJO_MIN)<=M%DY(J) .AND. &
            ABS( ((JJO_MAX-JJO_MIN+1)*MM%DY(JJO_MIN)-M%DY(J)) / MM%DY(JJO_MIN))>0.05 ) ALIGNED = .FALSE.
      IF ( (ABS(IOR)==1 .OR. ABS(IOR)==2) .AND. MM%DZ(KKO_MIN)<=M%DZ(K) .AND. &
            ABS( ((KKO_MAX-KKO_MIN+1)*MM%DZ(KKO_MIN)-M%DZ(K)) / MM%DZ(KKO_MIN))>0.05 ) ALIGNED = .FALSE.
      IF (.NOT.ALIGNED) THEN
         WRITE(LU_ERR,'(A,I3,A,I3)') 'ERROR: MESH ',NM,' is out of alignment with MESH ',NOM
         PROCESS_STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF

      SELECT CASE(ABS(IOR))
         CASE(1)
            IF ( (M%DY(J)<0.99_EB*MM%DY(JJO_MIN)) .AND. (M%DZ(K)>1.01_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
            IF ( (M%DY(J)>1.01_EB*MM%DY(JJO_MIN)) .AND. (M%DZ(K)<0.99_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
         CASE(2)
            IF ( (M%DX(I)<0.99_EB*MM%DX(IIO_MIN)) .AND. (M%DZ(K)>1.01_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
            IF ( (M%DX(I)>1.01_EB*MM%DX(IIO_MIN)) .AND. (M%DZ(K)<0.99_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
         CASE(3)
            IF ( (M%DY(J)<0.99_EB*MM%DY(JJO_MIN)) .AND. (M%DX(I)>1.01_EB*MM%DX(IIO_MIN)) ) ALIGNED = .FALSE.
            IF ( (M%DY(J)>1.01_EB*MM%DY(JJO_MIN)) .AND. (M%DX(I)<0.99_EB*MM%DX(IIO_MIN)) ) ALIGNED = .FALSE.
      END SELECT
      IF (.NOT.ALIGNED) THEN
         WRITE(LU_ERR,'(A,I3,A,I3)') 'ERROR: MESH ',NM,' is finer in one direction and coarser in the other than MESH ',NOM
         PROCESS_STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF

      M%IJKW(9,IW)  = NOM
      M%IJKW(10,IW) = IIO_MIN
      M%IJKW(11,IW) = JJO_MIN
      M%IJKW(12,IW) = KKO_MIN
      M%IJKW(13,IW) = IIO_MAX
      M%IJKW(14,IW) = JJO_MAX
      M%IJKW(15,IW) = KKO_MAX

      IF (I_OBST==0) THEN
         IF (.NOT.M%SOLID(ICG)) M%BOUNDARY_TYPE(IW) = INTERPOLATED_BOUNDARY
         IBCX = INTERPOLATED_SURF_INDEX
         M%IJKW(5,IW) = IBCX
      ENDIF
      IF (M%BOUNDARY_TYPE(IW) == INTERPOLATED_BOUNDARY) M%SOLID(M%CELL_INDEX(I,J,K)) = .FALSE.
      I_VENT = 0

   ENDIF FOUND_OTHER_MESH
 
ENDIF CHECK_MESHES

! Assign internal values of temp, density, and mass fraction
 
IF (N_SPECIES>0 .AND. M%BOUNDARY_TYPE(IW)/=VIRTUAL_BOUNDARY .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   M%RSUM(I,J,K) = M%RSUM(IIG,JJG,KKG)
   M%YY_F(IW,1:N_SPECIES)  = M%YY(IIG,JJG,KKG,1:N_SPECIES)
   M%YY(I,J,K,1:N_SPECIES) = M%YY(IIG,JJG,KKG,1:N_SPECIES)
ENDIF
 
M%RHO_F(IW) = M%RHO(IIG,JJG,KKG)
 
! Assign various other quantities to the cell
 
IF (I_OBST>0 .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   OBX=>M%OBSTRUCTION(I_OBST)
   M%AREA_ADJUST(IW) = OBX%INPUT_AREA(ABS(IOR))/OBX%FDS_AREA(ABS(IOR))
   IF (M%AREA_ADJUST(IW)==0._EB) M%AREA_ADJUST(IW) = 1._EB
   OBX%MASS = SURFACE(IBC)%SURFACE_DENSITY*M%AW(IW)*M%AREA_ADJUST(IW)
   IF (OBX%BULK_DENSITY>0._EB) OBX%MASS = OBX%BULK_DENSITY*(OBX%X2-OBX%X1)*(OBX%Y2-OBX%Y1)*(OBX%Z2-OBX%Z1)/OBX%VOLUME_ADJUST
ENDIF
 
! Prescribe exit velocity for surface cell
 
M%UW0(IW) = SURFACE(IBCX)%VEL

IF (I_OBST>0 .AND. SURFACE(IBCX)%VOLUME_FLUX/=0._EB .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   OBX=>M%OBSTRUCTION(I_OBST)
   M%UW0(IW) = SURFACE(IBCX)%VOLUME_FLUX/OBX%FDS_AREA(ABS(IOR))
ENDIF
IF (I_OBST>0 .AND. SURFACE(IBCX)%MASS_FLUX_TOTAL/=0._EB .AND. .NOT.EVACUATION_ONLY(NM)) THEN
   OBX=>M%OBSTRUCTION(I_OBST)
   M%UW0(IW) = SURFACE(IBCX)%MASS_FLUX_TOTAL / RHOA * M%AREA_ADJUST(IW)
ENDIF
 
T_ACTIVATE = T_BEGIN-1._EB
 
! Do VENT-specific set-ups
 
PROCESS_VENT: IF (I_VENT>0) THEN
 
   VT => M%VENTS(I_VENT)

   M%VENT_INDEX(IW) = I_VENT
   IF (.NOT.EVACUATION_ONLY(NM)) THEN
   M%AREA_ADJUST(IW) = VT%INPUT_AREA/VT%FDS_AREA
   IF (M%AREA_ADJUST(IW)==0._EB) M%AREA_ADJUST(IW) = 1._EB
   ENDIF
 
   ! Set the velocity at each surface cell
 
   M%UW0(IW) = SURFACE(IBCX)%VEL
 
   IF (SURFACE(IBCX)%VOLUME_FLUX    /=0._EB)   M%UW0(IW) = SURFACE(IBCX)%VOLUME_FLUX/VT%FDS_AREA
   IF (SURFACE(IBCX)%MASS_FLUX_TOTAL/=0._EB .AND. .NOT.EVACUATION_ONLY(NM)) &
        M%UW0(IW) = SURFACE(IBCX)%MASS_FLUX_TOTAL/RHOA*M%AREA_ADJUST(IW)
 
   ! Special velocity profiles
 
   IF (SURFACE(IBCX)%PROFILE==PARABOLIC) THEN 
      SELECT CASE(ABS(IOR))
         CASE(1)
            Y1 = M%Y(VT%J1)
            Y2 = M%Y(VT%J2)
            Z1 = M%Z(VT%K1)
            Z2 = M%Z(VT%K2)
            PY = 4._EB*(M%YC(J)-Y1)*(Y2-M%YC(J))/(Y2-Y1)**2
            PZ = 4._EB*(M%ZC(K)-Z1)*(Z2-M%ZC(K))/(Z2-Z1)**2
            M%UW0(IW) = M%UW0(IW)*PY*PZ
         CASE(2)
            X1 = M%X(VT%I1)
            X2 = M%X(VT%I2)
            Z1 = M%Z(VT%K1)
            Z2 = M%Z(VT%K2)
            PX = 4._EB*(M%XC(I)-X1)*(X2-M%XC(I))/(X2-X1)**2
            PZ = 4._EB*(M%ZC(K)-Z1)*(Z2-M%ZC(K))/(Z2-Z1)**2
            M%UW0(IW) = M%UW0(IW)*PX*PZ
         CASE(3)
            X1 = M%X(VT%I1)
            X2 = M%X(VT%I2)
            IF (CYLINDRICAL .AND. X1==0._EB) X1 = -X2
            Y1 = M%Y(VT%J1)
            Y2 = M%Y(VT%J2)
            PX = 4._EB*(M%XC(I)-X1)*(X2-M%XC(I))/(X2-X1)**2
            PY = 4._EB*(M%YC(J)-Y1)*(Y2-M%YC(J))/(Y2-Y1)**2
            IF (CYLINDRICAL) THEN
               M%UW0(IW) = M%UW0(IW)*PX
               ELSE
               M%UW0(IW) = M%UW0(IW)*PX*PY
               ENDIF
      END SELECT
   ENDIF
 
   IF (SURFACE(IBCX)%PROFILE==ATMOSPHERIC) THEN
      IF (M%ZC(K)<GROUND_LEVEL) THEN
         WRITE(LU_ERR,'(A,A,A)') 'ERROR: SURF ',TRIM(SURFACE(IBCX)%ID),' cannot be applied below GROUND_LEVEL'
         PROCESS_STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF
      M%UW0(IW) = M%UW0(IW)*((M%ZC(K)-GROUND_LEVEL)/SURFACE(IBCX)%Z0)**SURFACE(IBCX)%PLE
   ENDIF
 
   IF (SURFACE(IBCX)%PROFILE==ONED_PARABOLIC) THEN 
      SELECT CASE(ABS(IOR))
         CASE(1)
            Y1 = M%Y(VT%J1)
            Y2 = M%Y(VT%J2)
            PY = 4._EB*(M%YC(J)-Y1)*(Y2-M%YC(J))/(Y2-Y1)**2
            M%UW0(IW) = M%UW0(IW)*PY
         CASE(2)
            X1 = M%X(VT%I1)
            X2 = M%X(VT%I2)
            PX = 4._EB*(M%XC(I)-X1)*(X2-M%XC(I))/(X2-X1)**2
            M%UW0(IW) = M%UW0(IW)*PX
         CASE(3)
            X1 = M%X(VT%I1)
            X2 = M%X(VT%I2)
            Y1 = M%Y(VT%J1)
            Y2 = M%Y(VT%J2)
            PX = 4._EB*(M%XC(I)-X1)*(X2-M%XC(I))/(X2-X1)**2
            PY = 4._EB*(M%YC(J)-Y1)*(Y2-M%YC(J))/(Y2-Y1)**2
            M%UW0(IW) = M%UW0(IW)*PX*PY
      END SELECT
   ENDIF
 
   ! Check if fire spreads radially
 
   IF (VT%X0>-1.E5_EB) THEN
      DIST = SQRT((M%XW(IW)-VT%X0)**2 +(M%YW(IW)-VT%Y0)**2 +(M%ZW(IW)-VT%Z0)**2)
      T_ACTIVATE = T_BEGIN+DIST/VT%FIRE_SPREAD_RATE
   ENDIF

   ! Miscellaneous settings
 
   IF (VT%BOUNDARY_TYPE==OPEN_BOUNDARY   .AND. .NOT.M%SOLID(ICG)) M%BOUNDARY_TYPE(IW) = OPEN_BOUNDARY
   IF (VT%BOUNDARY_TYPE==MIRROR_BOUNDARY .AND. .NOT.M%SOLID(ICG)) M%BOUNDARY_TYPE(IW) = MIRROR_BOUNDARY

ENDIF PROCESS_VENT
 
! Set ignition time of each boundary cell
 
IF (T_ACTIVATE < T_BEGIN) THEN
   M%TW(IW) = SURFACE(IBCX)%T_IGN
ELSE
   M%TW(IW) = T_ACTIVATE
ENDIF
 
! Miscellaneous initializations

SF => SURFACE(IBCX)
M%E_WALL(IW) = SF%EMISSIVITY
IF (.NOT.EVACUATION_ONLY(NM)) M%NPPCW(IW)  = SF%NPPC      ! Number of particles per cell

! Initialize solid properties included in the WALL derived type variable
 
IF (M%BOUNDARY_TYPE(IW)/=NULL_BOUNDARY .OR. M%IJKW(9,IW)>0) CALL ALLOCATE_WALL_ARRAYS(NM,IW)

! Record original boundary condition index for exterior wall cells that might get covered up

IF (I_OBST==0) M%IBC_ORIG(IW) = IBCX

END SUBROUTINE INIT_WALL_CELL
 
 
SUBROUTINE OPEN_AND_CLOSE(T,NM)

! Check to see if a cell or OBSTruction is to be created or removed, or a VENT activated of deactivated

USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_STRINGS
USE CONTROL_VARIABLES, ONLY : CONTROL
USE DEVICE_VARIABLES, ONLY : DEVICE
REAL(EB), INTENT(IN) :: T
REAL(EB) :: T_BEGIN_TMP
INTEGER  :: N,II,JJ,KK,IW,IC
INTEGER, INTENT(IN) :: NM
LOGICAL :: CREATE_OBST,REMOVE_OBST,ACTIVATE_VENT,DEACTIVATE_VENT
CHARACTER(12) :: SV_LABEL
TYPE (VENTS_TYPE), POINTER :: VT

CALL POINT_TO_MESH(NM)
T_BEGIN_TMP = T_BEGIN
!IF (EVACUATION_ONLY(NM)) T_BEGIN_TMP = T_BEGIN - EVAC_DT_FLOWFIELD*EVAC_TIME_ITERATIONS
 
! Check to see if an obstacle is to be removed or created

OBST_LOOP: DO N=1,N_OBST
   OB=>OBSTRUCTION(N)
   IF (.NOT. OB%REMOVABLE) CYCLE OBST_LOOP
   CREATE_OBST = .FALSE.
   REMOVE_OBST = .FALSE.

   ! Over-ride DEVICE/CONTROL logic

   IF (OB%CONSUMABLE .AND. OB%MASS <= 0._EB) REMOVE_OBST = .TRUE.
   IF (OB%HIDDEN .AND. T<=T_BEGIN_TMP) THEN
      REMOVE_OBST = .TRUE.
      OB%HIDDEN = .FALSE.   ! The obstruction will not be removed if FDS thinks it is already HIDDEN
   ENDIF

   ! Decide if a DEVICE/CONTROL action is needed

   IF (.NOT.CREATE_OBST .AND. .NOT.REMOVE_OBST) THEN
      IF (OB%DEVC_INDEX > 0) THEN
         IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(OB%DEVC_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
         IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
            IF (.NOT.OB%HOLE_FILLER) CREATE_OBST = .TRUE.
            IF (     OB%HOLE_FILLER) REMOVE_OBST = .TRUE.
         ELSE
            IF (.NOT.OB%HOLE_FILLER) REMOVE_OBST = .TRUE.      
            IF (     OB%HOLE_FILLER) CREATE_OBST = .TRUE.      
         ENDIF
      ELSEIF (OB%CTRL_INDEX > 0) THEN
         IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(OB%CTRL_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
         IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
            IF (.NOT.OB%HOLE_FILLER) CREATE_OBST = .TRUE.
            IF (     OB%HOLE_FILLER) REMOVE_OBST = .TRUE.
         ELSE
            IF (.NOT.OB%HOLE_FILLER) REMOVE_OBST = .TRUE.      
            IF (     OB%HOLE_FILLER) CREATE_OBST = .TRUE.      
         ENDIF
      ENDIF
   ENDIF

  SV_LABEL  = 'null'

  IF (CREATE_OBST .AND. OB%HIDDEN) THEN
     OB%HIDDEN = .FALSE.
     SV_LABEL  = 'SHOW_OBST'
     CALL CREATE_OR_REMOVE_OBST(NM,OB%I1,OB%I2,OB%J1,OB%J2,OB%K1,OB%K2,1,N,T)
  ENDIF

   IF (REMOVE_OBST .AND. (.NOT. OB%HIDDEN)) THEN
      OB%HIDDEN = .TRUE.
      SV_LABEL  = 'HIDE_OBST'
      CALL CREATE_OR_REMOVE_OBST(NM,OB%I1,OB%I2,OB%J1,OB%J2,OB%K1,OB%K2,0,N,T)
   ENDIF

   ! Write a message to the Smokeview .smv file that the obstruction has been created or removed

   IF (SV_LABEL /= 'null') THEN
      IF (N_STRINGS+2>N_STRINGS_MAX) THEN
         CALL RE_ALLOCATE_STRINGS(NM)
         STRING => MESHES(NM)%STRING
      ENDIF
      N_STRINGS = N_STRINGS + 1
      WRITE(STRING(N_STRINGS),'(A,I3)') SV_LABEL,NM
      N_STRINGS = N_STRINGS + 1
      WRITE(STRING(N_STRINGS),'(I6,F10.2)') N,T
   ENDIF

ENDDO OBST_LOOP
 
! Check to see if a vent should be activated or deactivated
 
VENT_LOOP: DO N=1,N_VENT
   VT => VENTS(N)

   ACTIVATE_VENT   = .FALSE.
   DEACTIVATE_VENT = .FALSE.

   ! Over-ride DEVICE/CONTROL logic

   IF (.NOT.VT%ACTIVATED .AND. T<=T_BEGIN_TMP) DEACTIVATE_VENT = .TRUE.

   ! Decide if a VENT is to activate or de-activate based on a DEVICE or CONTROLLER

   IF (.NOT.ACTIVATE_VENT .AND. .NOT.DEACTIVATE_VENT) THEN
      IF (VT%DEVC_INDEX > 0) THEN
         IF (DEVICE(VT%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(VT%DEVC_INDEX)%PRIOR_STATE) CYCLE VENT_LOOP
         IF (DEVICE(VT%DEVC_INDEX)%CURRENT_STATE) THEN
            ACTIVATE_VENT   = .TRUE.
         ELSE
            DEACTIVATE_VENT = .TRUE.
         ENDIF
      ELSEIF (VT%CTRL_INDEX > 0) THEN
         IF (CONTROL(VT%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(VT%CTRL_INDEX)%PRIOR_STATE) CYCLE VENT_LOOP
         IF (CONTROL(VT%CTRL_INDEX)%CURRENT_STATE) THEN
            ACTIVATE_VENT   = .TRUE.
         ELSE
            DEACTIVATE_VENT = .TRUE.
         ENDIF
      ENDIF
   ENDIF

   IF (.NOT.ACTIVATE_VENT .AND. .NOT.DEACTIVATE_VENT) CYCLE VENT_LOOP

   ! Find the wall indices (IW) for the vent and set the activation time (TW) for each one

   DO KK=VT%K1+1,MAX(VT%K1+1,VT%K2)
      DO JJ=VT%J1+1,MAX(VT%J1+1,VT%J2)
         DO II=VT%I1+1,MAX(VT%I1+1,VT%I2)
            SELECT CASE(VT%IOR)
               CASE(1:)
                  IC = CELL_INDEX(II,JJ,KK)
               CASE(-1)
                  IC = CELL_INDEX(II-1,JJ,KK)
               CASE(-2)
                  IC = CELL_INDEX(II,JJ-1,KK)
               CASE(-3)
                  IC = CELL_INDEX(II,JJ,KK-1)
            END SELECT
            IW = WALL_INDEX(IC,-VT%IOR)
            IF (IW==0) CYCLE

            IF (ACTIVATE_VENT) THEN
               IF (VT%X0>-1.E5_EB) THEN
                  TW(IW) = T + SQRT((XW(IW)-VT%X0)**2 +(YW(IW)-VT%Y0)**2 +(ZW(IW)-VT%Z0)**2)/VT%FIRE_SPREAD_RATE
               ELSE
                  TW(IW) = T        
               ENDIF
            ELSE
               TW(IW) = 1000000._EB
            ENDIF
         ENDDO
      ENDDO
   ENDDO

   ! Write message to .smv file

   IF (ACTIVATE_VENT)   SV_LABEL = 'OPEN_VENT'
   IF (DEACTIVATE_VENT) SV_LABEL = 'CLOSE_VENT'

   IF (N_STRINGS+2>N_STRINGS_MAX) THEN
      CALL RE_ALLOCATE_STRINGS(NM)
      STRING => MESHES(NM)%STRING
   ENDIF
   N_STRINGS = N_STRINGS + 1
   WRITE(STRING(N_STRINGS),'(A,I3)') SV_LABEL,NM
   N_STRINGS = N_STRINGS + 1
   WRITE(STRING(N_STRINGS),'(I6,F10.2)') N,T

ENDDO VENT_LOOP
 
END SUBROUTINE OPEN_AND_CLOSE
 

 
SUBROUTINE CREATE_OR_REMOVE_OBST(NM,I1,I2,J1,J2,K1,K2,CR_INDEX,OBST_INDEX,T)

! Create or remove the obstruction whose NODES (not cells) are given by I1, I2, etc.

USE GEOMETRY_FUNCTIONS, ONLY : BLOCK_CELL
INTEGER :: I1,I2,J1,J2,K1,K2,I,J,K,IW,ICG,IC
INTEGER, INTENT(IN) :: NM,CR_INDEX,OBST_INDEX
REAL(EB) :: T
LOGICAL :: CREATE,REMOVE
 
CALL POINT_TO_MESH(NM)

REMOVE = .FALSE.
CREATE = .FALSE.
IF (CR_INDEX==0) REMOVE = .TRUE.
IF (CR_INDEX==1) CREATE = .TRUE.
 
! Blank or unblank cells that make up the OBSTruction
 
CALL BLOCK_CELL(NM,I1+1,I2,J1+1,J2,K1+1,K2,CR_INDEX,OBST_INDEX)
 
! Process the x boundaries of the OBSTruction

DO K=K1+1,K2
   DO J=J1+1,J2
      IC  = CELL_INDEX(I1+1,J,K)
      ICG = CELL_INDEX(I1  ,J,K)
      IW  = WALL_INDEX(ICG,1)
      IF (I1>0)    CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I1  ,J,K)
      ICG = CELL_INDEX(I1+1,J,K)
      IW  = WALL_INDEX(ICG,-1)
      IF (I1<IBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I2  ,J,K)
      ICG = CELL_INDEX(I2+1,J,K)
      IW  = WALL_INDEX(ICG,-1)
      IF (I2<IBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I2+1,J,K)
      ICG = CELL_INDEX(I2  ,J,K)
      IW  = WALL_INDEX(ICG, 1)
      IF (I2>0)    CALL GET_BOUNDARY_TYPE
   ENDDO
ENDDO

! Process the y boundaries of the OBSTruction

DO K=K1+1,K2
   DO I=I1+1,I2
      IC  = CELL_INDEX(I,J1+1,K)
      ICG = CELL_INDEX(I,J1  ,K)
      IW  = WALL_INDEX(ICG, 2)
      IF (J1>0)    CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J1  ,K)
      ICG = CELL_INDEX(I,J1+1,K)
      IW  = WALL_INDEX(ICG,-2)
      IF (J1<JBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J2  ,K)
      ICG = CELL_INDEX(I,J2+1,K)
      IW  = WALL_INDEX(ICG,-2)
      IF (J2<JBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J2+1,K)
      ICG = CELL_INDEX(I,J2  ,K)
      IW  = WALL_INDEX(ICG, 2)
      IF (J2>0)    CALL GET_BOUNDARY_TYPE
   ENDDO
ENDDO

! Process the z boundaries of the OBSTruction

DO J=J1+1,J2
   DO I=I1+1,I2
      IC  = CELL_INDEX(I,J,K1+1)
      ICG = CELL_INDEX(I,J,K1  )
      IW  = WALL_INDEX(ICG, 3)
      IF (K1>0)    CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J,K1  )
      ICG = CELL_INDEX(I,J,K1+1)
      IW  = WALL_INDEX(ICG,-3)
      IF (K1<KBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J,K2  )
      ICG = CELL_INDEX(I,J,K2+1)
      IW  = WALL_INDEX(ICG,-3)
      IF (K2<KBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J,K2+1)
      ICG = CELL_INDEX(I,J,K2  )
      IW  = WALL_INDEX(ICG, 3)
      IF (K2>0)    CALL GET_BOUNDARY_TYPE
   ENDDO
ENDDO


! Nullify block edges on blockage that is to be removed
 
DO K=K1,K2
   DO J=J1,J2
      IF (J>J1) CALL REDEFINE_EDGE(I1,J,K,2)
      IF (J>J1) CALL REDEFINE_EDGE(I2,J,K,2)
      IF (K>K1) CALL REDEFINE_EDGE(I1,J,K,3)
      IF (K>K1) CALL REDEFINE_EDGE(I2,J,K,3)
   ENDDO
ENDDO

DO K=K1,K2
   DO I=I1,I2
      IF (I>I1) CALL REDEFINE_EDGE(I,J1,K,1)
      IF (I>I1) CALL REDEFINE_EDGE(I,J2,K,1)
      IF (K>K1) CALL REDEFINE_EDGE(I,J1,K,3)
      IF (K>K1) CALL REDEFINE_EDGE(I,J2,K,3)
   ENDDO
ENDDO

DO J=J1,J2
   DO I=I1,I2
      IF (I>I1) CALL REDEFINE_EDGE(I,J,K1,1)
      IF (I>I1) CALL REDEFINE_EDGE(I,J,K2,1)
      IF (J>J1) CALL REDEFINE_EDGE(I,J,K1,2)
      IF (J>J1) CALL REDEFINE_EDGE(I,J,K2,2)
   ENDDO
ENDDO

CONTAINS

SUBROUTINE GET_BOUNDARY_TYPE

USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_WALL_ARRAYS
INTEGER :: I_OBST,IBC,IOR,IIG,JJG,KKG,IW_OLD,IBC_OLD

IF (IW==0) RETURN

IBC = IJKW(5,IW)

IF (IW<=NEWC .AND. REMOVE) THEN
   BOUNDARY_TYPE(IW) = SOLID_BOUNDARY
   IF (IBC_ORIG(IW)==MIRROR_SURF_INDEX)       BOUNDARY_TYPE(IW) = MIRROR_BOUNDARY
   IF (IBC_ORIG(IW)==OPEN_SURF_INDEX)         BOUNDARY_TYPE(IW) = OPEN_BOUNDARY
   IF (IBC_ORIG(IW)==INTERPOLATED_SURF_INDEX) THEN
      BOUNDARY_TYPE(IW) = INTERPOLATED_BOUNDARY
   ENDIF
ENDIF

IF (IW>NEWC) THEN
   I_OBST = OBST_INDEX_W(IW)
   IF (I_OBST>0 .AND.      OBSTRUCTION(I_OBST)%HIDDEN .AND. .NOT.SOLID(IC) ) BOUNDARY_TYPE(IW) = NULL_BOUNDARY
   IF (I_OBST>0 .AND. .NOT.OBSTRUCTION(I_OBST)%HIDDEN .AND. .NOT.SOLID(ICG)) BOUNDARY_TYPE(IW) = SOLID_BOUNDARY
   IF (     SOLID(ICG)) BOUNDARY_TYPE(IW) = NULL_BOUNDARY
ENDIF

IF (CREATE) THEN
   IF (SOLID(ICG)) THEN
      BOUNDARY_TYPE(IW) = NULL_BOUNDARY
   ELSE
      BOUNDARY_TYPE(IW) = SOLID_BOUNDARY
      IF (SURFACE(IBC)%POROUS) BOUNDARY_TYPE(IW) = POROUS_BOUNDARY
      IF (TW(IW)<T) TW(IW)=T
   ENDIF
ENDIF

! Make the ghost cell either solid or not depending on the boundary type

IF (IW<=NEWC .AND. (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY .OR. BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY)) SOLID(IC) = .FALSE.
IF (IW<=NEWC .AND. (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==SOLID_BOUNDARY))        SOLID(IC) = .TRUE.

! Initialize solid properties included in the WALL derived type variable
 
IF (BOUNDARY_TYPE(IW)/=NULL_BOUNDARY) CALL ALLOCATE_WALL_ARRAYS(NM,IW)

! Inherit ignition time from BURNed_AWAY obstruction

IF (REMOVE .AND. SURFACE(IBC)%PYROLYSIS_MODEL==PYROLYSIS_SPECIFIED) THEN
   IIG = IJKW(6,IW)
   JJG = IJKW(7,IW)
   KKG = IJKW(8,IW)
   IOR = IJKW(4,IW)
   IW_OLD = 0
   IC     = 0
   SELECT CASE(IOR)
      CASE(-1)
         IF (IIG>1)    IC = CELL_INDEX(IIG-1,JJG,KKG)
      CASE( 1)
         IF (IIG<IBAR) IC = CELL_INDEX(IIG+1,JJG,KKG)
      CASE(-2)
         IF (JJG>1)    IC = CELL_INDEX(IIG,JJG-1,KKG)
      CASE( 2)
         IF (JJG<JBAR) IC = CELL_INDEX(IIG,JJG+1,KKG)
      CASE(-3)
         IF (KKG>1)    IC = CELL_INDEX(IIG,JJG,KKG-1)
      CASE( 3)
         IF (KKG<KBAR) IC = CELL_INDEX(IIG,JJG,KKG+1)
   END SELECT
   IW_OLD = WALL_INDEX(IC,-IOR)
   IF (IW_OLD>0) THEN
      IBC_OLD = IJKW(5,IW_OLD)
      IF (IBC==IBC_OLD) TW(IW) = TW(IW_OLD)
   ENDIF
ENDIF

END SUBROUTINE GET_BOUNDARY_TYPE

SUBROUTINE REDEFINE_EDGE(II,JJ,KK,IEC)

! Change a few properties of the EDGEs that have been exposed or covered up by a blockage

INTEGER :: IE,II,JJ,KK,IEC

SELECT CASE(IEC)
   CASE(1)
      IE = EDGE_INDEX(CELL_INDEX(II,JJ,KK), 4)
   CASE(2)
      IE = EDGE_INDEX(CELL_INDEX(II,JJ,KK), 8)
   CASE(3)
      IE = EDGE_INDEX(CELL_INDEX(II,JJ,KK),12)
END SELECT

IF (IE>0) THEN
   IF (REMOVE .AND. IJKE( 9,IE)/=0) EDGE_TYPE(IE,1) = INTERPOLATED_EDGE
   IF (REMOVE .AND. IJKE(13,IE)/=0) EDGE_TYPE(IE,2) = INTERPOLATED_EDGE
   IF (CREATE .AND. IJKE( 9,IE)/=0) EDGE_TYPE(IE,1) = SOLID_EDGE
   IF (CREATE .AND. IJKE(13,IE)/=0) EDGE_TYPE(IE,2) = SOLID_EDGE
ENDIF

END SUBROUTINE REDEFINE_EDGE

END SUBROUTINE CREATE_OR_REMOVE_OBST


SUBROUTINE DEFINE_EDGE(II,JJ,KK,IOR,IEC,NM,I_OBST,IERR)
 
! Set up edge arrays for velocity boundary conditions
 
INTEGER, INTENT(IN) :: II,JJ,KK,IOR,IEC,NM
INTEGER :: NOM,ICMM,ICMP,ICPM,ICPP,I_OBST,IE,IW,IIO,JJO,KKO,IW1,IW2,IERR
REAL(EB) :: XI,YJ,ZK
TYPE (MESH_TYPE), POINTER :: MM
 
IF (I_OBST>0) OB=>OBSTRUCTION(I_OBST)
 
! Find the wall cells on each side of the edge 

IW1 = -1
IW2 = -1
 
EDGE_DIRECTION_1: SELECT CASE(IEC)
   CASE(1) EDGE_DIRECTION_1
      SELECT CASE(IOR)
         CASE(-2)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ,KK)  ,2)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ,KK+1),2)
         CASE( 2)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK)  ,-2)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK+1),-2)
         CASE(-3)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ  ,KK),3)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK),3)
         CASE( 3)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ  ,KK+1),-3)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK+1),-3)
      END SELECT
   CASE(2) EDGE_DIRECTION_1
      SELECT CASE(IOR)
         CASE(-1)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ,KK)  ,1)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ,KK+1),1)
         CASE( 1)
            IW1 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK)  ,-1)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK+1),-1)
         CASE(-3)
            IW1 = WALL_INDEX(CELL_INDEX(II  ,JJ,KK),3)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK),3)
         CASE( 3)
            IW1 = WALL_INDEX(CELL_INDEX(II  ,JJ,KK+1),-3)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK+1),-3)
      END SELECT
   CASE(3) EDGE_DIRECTION_1
      SELECT CASE(IOR)
         CASE(-1)
            IW1 = WALL_INDEX(CELL_INDEX(II,JJ  ,KK),1)
            IW2 = WALL_INDEX(CELL_INDEX(II,JJ+1,KK),1)
         CASE( 1)
            IW1 = WALL_INDEX(CELL_INDEX(II+1,JJ  ,KK),-1)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ+1,KK),-1)
         CASE(-2)
            IW1 = WALL_INDEX(CELL_INDEX(II  ,JJ,KK),2)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ,KK),2)
         CASE( 2)
            IW1 = WALL_INDEX(CELL_INDEX(II  ,JJ+1,KK),-2)
            IW2 = WALL_INDEX(CELL_INDEX(II+1,JJ+1,KK),-2)
      END SELECT
END SELECT EDGE_DIRECTION_1

! Decide what to do based on whether or not adjacent tiles exist

IF (IW1==-1 .OR. IW2==-1) THEN
   WRITE(LU_ERR,'(A,I2,A,3I3)') 'ERROR: Edge initialization failed; Mesh: ',NM,', Cell: ',II,JJ,KK
   PROCESS_STOP_STATUS = SETUP_STOP
   IERR = 1
   RETURN
ENDIF
IF (IW1==0 .AND. IW2==0) RETURN
IF (IW1> 0 .AND. IW2==0) IW = IW1
IF (IW1==0 .AND. IW2> 0) IW = IW2
IF (IW1> 0 .AND. IW2> 0) IW = IW2  ! Arbitrary decision
 
! Assign the Index of the Edge (IE) and add to the list

ICMM = CELL_INDEX(II,JJ,KK)
SELECT CASE(IEC)
   CASE(1)
      IE = EDGE_INDEX(ICMM,4)
   CASE(2)
      IE = EDGE_INDEX(ICMM,8)
   CASE(3)
      IE = EDGE_INDEX(ICMM,12)
END SELECT
 
IF (IE==0) THEN
   N_EDGES = N_EDGES + 1
   IE = N_EDGES
ENDIF
 
! Determine the wall index of the adjacent wall tile
 
NOM = 0
IIO = 0
JJO = 0
KKO = 0
 
IF (IJKW(9,IW)>0) THEN
   NOM = IJKW( 9,IW)
   IIO = IJKW(10,IW)
   JJO = IJKW(11,IW)
   KKO = IJKW(12,IW)
ENDIF

! Fill up array IJKE with edge parameters
 
IJKE(1,IE) = II
IJKE(2,IE) = JJ
IJKE(3,IE) = KK
IJKE(4,IE) = IEC
 
! Fill in EDGE_INDEX and the rest of IJKE

EDGE_DIRECTION_2: SELECT CASE(IEC)
 
   CASE (1) EDGE_DIRECTION_2
 
      IF (ABS(IOR)==2 .AND. BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) EDGE_TYPE(IE,1) = INTERPOLATED_EDGE
      IF (ABS(IOR)==3 .AND. BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) EDGE_TYPE(IE,2) = INTERPOLATED_EDGE
      ICPM = CELL_INDEX(II,JJ+1,KK)
      ICPP = CELL_INDEX(II,JJ+1,KK+1)
      ICMP = CELL_INDEX(II,JJ,KK+1)
      IJKE(5,IE) = ICMM
      IJKE(6,IE) = ICPM
      IJKE(7,IE) = ICMP
      IJKE(8,IE) = ICPP
      EDGE_INDEX(ICPP,1) = IE
      EDGE_INDEX(ICMP,2) = IE
      EDGE_INDEX(ICPM,3) = IE
      EDGE_INDEX(ICMM,4) = IE
      IF (NOM/=0) THEN
         SELECT CASE(ABS(IOR))
            CASE(2)
               IF (IOR>0) IJKE( 9,IE) = -NOM
               IF (IOR<0) IJKE( 9,IE) =  NOM
               IJKE(10,IE) = IIO
               IJKE(11,IE) = JJO
               MM => MESHES(NOM)
               ZK  = MIN( REAL(MM%KBAR,EB)+ALMOST_ONE , MM%CELLSK(NINT((Z(KK)-MM%ZS)*MM%RDZINT))+1._EB )
               KKO = MAX(1,FLOOR(ZK))
               M%EDGE_INTERPOLATION_FACTOR(IE,1) = ZK-KKO
               IJKE(12,IE) = KKO
   
            CASE(3)
               IF (IOR>0) IJKE(13,IE) = -NOM
               IF (IOR<0) IJKE(13,IE) =  NOM
               IJKE(14,IE) = IIO
               MM => MESHES(NOM)
               YJ  = MIN( REAL(MM%JBAR,EB)+ALMOST_ONE , MM%CELLSJ(NINT((Y(JJ)-MM%YS)*MM%RDYINT))+1._EB )
               JJO = MAX(1,FLOOR(YJ))
               M%EDGE_INTERPOLATION_FACTOR(IE,2) = YJ-JJO
               IJKE(15,IE) = JJO
               IJKE(16,IE) = KKO
         END SELECT
      ENDIF
 
   CASE (2) EDGE_DIRECTION_2
    
      IF (ABS(IOR)==3 .AND. BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) EDGE_TYPE(IE,1) = INTERPOLATED_EDGE
      IF (ABS(IOR)==1 .AND. BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) EDGE_TYPE(IE,2) = INTERPOLATED_EDGE
      ICMP = CELL_INDEX(II+1,JJ,KK)
      ICPP = CELL_INDEX(II+1,JJ,KK+1)
      ICPM = CELL_INDEX(II,JJ,KK+1)
      IJKE(5,IE) = ICMM
      IJKE(6,IE) = ICPM
      IJKE(7,IE) = ICMP
      IJKE(8,IE) = ICPP
      EDGE_INDEX(ICPP,5) = IE
      EDGE_INDEX(ICPM,6) = IE
      EDGE_INDEX(ICMP,7) = IE
      EDGE_INDEX(ICMM,8) = IE
      IF (NOM/=0) THEN
         SELECT CASE(ABS(IOR))
            CASE( 1)
               IF (IOR>0) IJKE(13,IE) = -NOM
               IF (IOR<0) IJKE(13,IE) =  NOM
               IJKE(14,IE) = IIO
               IJKE(15,IE) = JJO
               MM => MESHES(NOM)
               ZK  = MIN( REAL(MM%KBAR,EB)+ALMOST_ONE , MM%CELLSK(NINT((Z(KK)-MM%ZS)*MM%RDZINT))+1._EB )
               KKO = MAX(1,FLOOR(ZK))
               M%EDGE_INTERPOLATION_FACTOR(IE,2) = ZK-KKO
               IJKE(16,IE) = KKO
            CASE( 3)
               IF (IOR>0) IJKE( 9,IE) = -NOM
               IF (IOR<0) IJKE( 9,IE) =  NOM
               MM => MESHES(NOM)
               XI  = MIN( REAL(MM%IBAR,EB)+ALMOST_ONE , MM%CELLSI(NINT((X(II)-MM%XS)*MM%RDXINT))+1._EB )
               IIO = MAX(1,FLOOR(XI))
               M%EDGE_INTERPOLATION_FACTOR(IE,1) = XI-IIO
               IJKE(10,IE) = IIO
               IJKE(11,IE) = JJO
               IJKE(12,IE) = KKO
         END SELECT 
      ENDIF
    
   CASE (3) EDGE_DIRECTION_2
 
      IF (ABS(IOR)==1 .AND. BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) EDGE_TYPE(IE,1) = INTERPOLATED_EDGE
      IF (ABS(IOR)==2 .AND. BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) EDGE_TYPE(IE,2) = INTERPOLATED_EDGE
      ICPM = CELL_INDEX(II+1,JJ,KK)
      ICPP = CELL_INDEX(II+1,JJ+1,KK)
      ICMP = CELL_INDEX(II,JJ+1,KK)
      IJKE(5,IE) = ICMM
      IJKE(6,IE) = ICPM
      IJKE(7,IE) = ICMP
      IJKE(8,IE) = ICPP
      EDGE_INDEX(ICPP, 9) = IE
      EDGE_INDEX(ICMP,10) = IE
      EDGE_INDEX(ICPM,11) = IE
      EDGE_INDEX(ICMM,12) = IE
      IF (NOM/=0) THEN
         SELECT CASE(ABS(IOR))
            CASE( 1)
               IF (IOR>0) IJKE( 9,IE) = -NOM
               IF (IOR<0) IJKE( 9,IE) =  NOM
               IJKE(10,IE) = IIO
               MM => MESHES(NOM)
               YJ  = MIN( REAL(MM%JBAR,EB)+ALMOST_ONE , MM%CELLSJ(NINT((Y(JJ)-MM%YS)*MM%RDYINT))+1._EB )
               JJO = MAX(1,FLOOR(YJ))
               M%EDGE_INTERPOLATION_FACTOR(IE,1) = YJ-JJO
               IJKE(11,IE) = JJO
               IJKE(12,IE) = KKO
            CASE( 2)
               IF (IOR>0) IJKE(13,IE) = -NOM
               IF (IOR<0) IJKE(13,IE) =  NOM
               MM => MESHES(NOM)
               XI  = MIN( REAL(MM%IBAR,EB)+ALMOST_ONE , MM%CELLSI(NINT((X(II)-MM%XS)*MM%RDXINT))+1._EB )
               IIO = MAX(1,FLOOR(XI))
               M%EDGE_INTERPOLATION_FACTOR(IE,2) = XI-IIO
               IJKE(14,IE) = IIO
               IJKE(15,IE) = JJO
               IJKE(16,IE) = KKO
         END SELECT
      ENDIF
 
END SELECT EDGE_DIRECTION_2

! Special "free-slip" cases

IF (I_OBST>0) THEN
   IF (.NOT.OB%SAWTOOTH) EDGE_TYPE(IE,:) = SMOOTH_EDGE
ENDIF

IF (EVACUATION_ONLY(NM)) EDGE_TYPE(IE,:) = SMOOTH_EDGE

END SUBROUTINE DEFINE_EDGE



SUBROUTINE INITIAL_NOISE(NM)

! Generate random noise at the start of the simulation
 
REAL(EB) :: VFAC,RN
INTEGER  :: I,J,K
INTEGER, INTENT(IN) :: NM
 
IF (EVACUATION_ONLY(NM)) RETURN
 
CALL POINT_TO_MESH(NM)

VFAC = 0.005_EB
 
DO K=0,KBAR
   DO J=0,JBAR
      DO I=1,IBAR
         IF (SOLID(CELL_INDEX(I,J,K))   .OR. SOLID(CELL_INDEX(I,J,K+1)) .OR. &
             SOLID(CELL_INDEX(I,J+1,K)) .OR. SOLID(CELL_INDEX(I,J+1,K+1)))  CYCLE 
         CALL RANDOM_NUMBER(RN)
         RN = VFAC*(-1._EB + 2._EB*RN)*CELL_SIZE
         W(I,J,K)   = W(I,J,K)   - RN*RDY(J)
         W(I,J+1,K) = W(I,J+1,K) + RN*RDY(J+1)
         V(I,J,K)   = V(I,J,K)   + RN*RDZ(K)
         V(I,J,K+1) = V(I,J,K+1) - RN*RDZ(K+1)
      ENDDO
   ENDDO
ENDDO
DO K=0,KBAR
   DO J=1,JBAR
      DO I=0,IBAR
         IF (SOLID(CELL_INDEX(I,J,K))   .OR. SOLID(CELL_INDEX(I,J,K+1)) .OR. &
             SOLID(CELL_INDEX(I+1,J,K)) .OR. SOLID(CELL_INDEX(I+1,J,K+1)))  CYCLE 
         CALL RANDOM_NUMBER(RN)
         RN = VFAC*(-1._EB + 2._EB*RN)*CELL_SIZE
         W(I,J,K)   = W(I,J,K)   - RN*RDX(I)*R(I)*RRN(I)
         W(I+1,J,K) = W(I+1,J,K) + RN*RDX(I+1)*R(I)*RRN(I+1)
         U(I,J,K)   = U(I,J,K)   + RN*RDZ(K)
         U(I,J,K+1) = U(I,J,K+1) - RN*RDZ(K+1)
      ENDDO 
   ENDDO
ENDDO
DO K=1,KBAR
   DO J=0,JBAR
      DO I=0,IBAR
         IF (SOLID(CELL_INDEX(I,J,K))   .OR. SOLID(CELL_INDEX(I,J+1,K)) .OR. &
             SOLID(CELL_INDEX(I+1,J,K)) .OR. SOLID(CELL_INDEX(I+1,J+1,K)))  CYCLE
         CALL RANDOM_NUMBER(RN)
         RN = VFAC*(-1._EB + 2._EB*RN)*CELL_SIZE
         V(I,J,K)   = V(I,J,K)   - RN*RDX(I)
         V(I+1,J,K) = V(I+1,J,K) + RN*RDX(I+1)
         U(I,J,K)   = U(I,J,K)   + RN*RDY(J)
         U(I,J+1,K) = U(I,J+1,K) - RN*RDY(J+1)
      ENDDO 
   ENDDO
ENDDO

END SUBROUTINE INITIAL_NOISE
 

SUBROUTINE GET_REV_init(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE

WRITE(MODULE_DATE,'(A)') initrev(INDEX(initrev,':')+1:LEN_TRIM(initrev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') initdate

END SUBROUTINE GET_REV_init

 
END MODULE INIT
